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ABSTRACT 


In this research, a novel approach for conducting signals intelligence from a 
distributed network of wireless nodes is developed. The primary objective of this 
research is enhancing signal collection in a specified target direction. Two conflicting 
priorities are addressed. One is the time required to determine the target direction and 
form the beams. The other is the energy consumption involved in developing these 
solutions. 

Two competing enhanced collection methodologies (ECM), ECM-1 and ECM-2, 
were developed and analyzed. ECM-1 uses a combination of time difference of arrival 
(TDOA) and adaptive beamforming. ECM-2 uses adaptive beamforming that performs a 
beamscan similar to phased-array radars. Additionally, two competing methods for 
forming the beams are developed. Method One uses data exclusively from the same 
elements. Method Two uses data from a new subset of sensors, for each iteration. 

Analytical expressions were derived for energy consumption and the time 
required to develop, to compare the competing methodologies. ECM-1 is shown to be far 
superior to ECM-2 in both energy consumption and the time required to enhance signal 
collection, whereas Method Two is shown to be far superior to Method One in formation 
of the beams. 
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EXECUTIVE SUMMARY 


In this research, a novel approach for conducting signals intelligence from a 
distributed network of wireless nodes was developed. The primary objective of this 
research is enhancing collection in a specified target direction. The SIGINT/IW wireless 
sensor network consists of a large number of small sensor nodes that are distributed over 
a large operational area, to acquire information about signals of interest. These sensor 
nodes collaborate among themselves to form an ad-hoc network and disseminate the 
collected target information to a centralized controller node, e.g., UAV, UUV, etc. In this 
research, it is assumed that the clustering functions to create a sensor grid have been 
completed. It is assumed that perfect node location has been determined and that these 
data are available to the centralized controller node. 

The centralized controller node, henceforth referred to as the central controller, 
acts as a reference antenna by gaining initial intercept and frequency of interest 
determination. The central controller is also tasked with maintaining frequency, phase 
and data synchronization among the nodes. Since it is envisioned that the sensors within 
this network will utilize low cost, omnidirectional antennas, adaptive beamforming is 
used in order to enhance the RF signal collection and/or information operations 
capability. 

The task of conducting RF collection from a network of this type was divided into 
four steps. The first step is the determination of the frequency of interest of the target. 

The second step is the synchronization of the individual nodes. The third step is the 
determination of the direction of the target. Finally, the fourth step is the formation of the 
beam in the target direction. Two of these steps are specifically addressed in this 
research, steps three and four. The first step will be assumed to be accomplished by an 
implementation of the periodogram. The second step will be assumed to be accomplished 
through a “brute-force” method similar to the one presented by, Jenn, et al. [19] for 
wireless radar applications. 
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When addressing steps three and four, two eonflieting priorities were addressed. 
One is the time required to determine the target direetion and the time required to form 
the beams. The other is the energy consumption involved in developing these solutions. 
Two competing enhanced collection methodologies (ECM), called ECM-1 and ECM-2, 
were developed and researched. ECM-1 uses a combination of time difference of arrival 
(TDOA) and adaptive beamforming. ECM-2 uses adaptive beamforming that performs a 
beamscan similar to phased-array radars. 

Additionally, two competing methods for forming the beams were created. 

Method One and Method Two. Method One created a subset of sensors from a randomly 
distributed antenna array. This method uses data exclusively from these elements and 
increases the number of iterations until the desired beam is formed. Method Two creates 
a new subset of sensors, from a randomly distributed array, for each iteration, until the 
desired beam is formed. 

Analytical expressions were derived that proved that Method Two provides better 
power management across the sensor network. Additionally, Method Two was able to 
mitigate the grating lobes that are present in arrays where the antenna element spacing is 
greater than a half-wavelength. 

ECM-1 is described as follows. The central controller will utilize the cross¬ 
ambiguity function to calculate the TDOAs between the individual sensors and the 
central controller. The central controller will then use the Newton-Raphson technique 
(NRT) to determine a line of bearing to the target emitter. The central controller using 
Method Two will begin the beamforming process. The central controller randomly selects 
the number of elements required to create the appropriate beam for a given azimuth and 
elevation of a target emitter. 

ECM-2 is described as follows. The central controller will then coordinate the 
direction-of-arrival (DOA) determination using beamscanning. Eor this method, the 
central controller employs Method Two for creating the beams. The central controller 
randomly selects the number of array elements needed to form the desired beam. This 
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process is repeated until a 360 degree scan has been completed. The central controller 
will then make a determination as to the direetion of the highest reeeived power. 

Analytical expressions, for energy consumption and time to develop a solution, 
were derived to compare ECM-1 to ECM-2. ECM-1 was shown to be far superior to 
ECM-2 in both energy consumption and the time required to formulate a solution. 

In this research, a novel idea for conducting signals collection from a wireless 
sensor network was developed. To facilitate signals collection from this network, an 
enhanced collection methodology, ECM-1, was developed. ECM-1 is comprised of a 
unique rapid, energy efficient technique for determining the direction of the target and a 
unique adaptive beamforming technique. Method Two. Analytical expressions were 
derived to compare ECM-1 with ECM-2. Several simulations were developed and 
implemented to study the differences between ECM-1 and ECM-2. 
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I. INTRODUCTION 


In this research, a novel idea for conducting signals collection from a wireless 
sensor network was developed. Our current national intelligence resources are limited in 
responsiveness and number. Current tactical intelligence collection sensors lack 
capability, flexibility and are limited by the inability to connect with national sensors. 
Both Joint Vision 2020 and Sea Power recognize the role of intelligence collection 
networks in full-spectrum dominance and the need for persistent intelligence, 
surveillance and reconnaissance operations [20, 23]. Advances in software radio, wireless 
communications and electronics have enabled the development of smart, disposable 
microsensors that can be deployed in the battlespace for conducting Signals Intelligence 
(SIGINT) and Information Warfare (IW) [2, 73]. These improvements in sensor network 
technologies continue to decrease the size, weight and costs of sensors and increase their 
resolution and utility [9, 65]. 

Previously, wireless sensors focused on the collection of data with low bandwidth 
requirements. The methodology in this research will require and assume that sufficient 
bandwidth is available to transport the required information back to a centralized node for 
processing. This research will leverage and expand on current research areas: 
“opportunistic array” research by Jenn, et al. [20], a second is adaptive beamforming for 
sensor networks by Tummala, McEachen, Vincent, et al. [68], and another is time 
difference of arrival estimation by Loomis, et al. [29]. 

The “opportunistic array” concept is an integrated digital phased array radar, 
where the array elements connect wirelessly and are placed randomly over the surface of 
a ship. The unique requirement of having these elements function coherently led Jenn, et 
al. [21]. to design individual transmit and receive (T/R) modules. Specifically, these 
researchers focused on the wireless local oscillator (LO) distribution and control 
methods, microstrip patch antenna design, commercial-off-the-shelf (COTS) hardware 
investigation and synchronization. 

In a separate research project at the Naval Postgraduate School, Vincent, et al 
explored the concept of beamforming in order to minimize power requirements for data 
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exfiltration. The results were positive analytieally proving that by treating a sensor grid as 
an array of randomly distributed antenna elements, a beamforming algorithm eould be 
used to enhanee the transmission range of the network. 

Building and testing a rapidly deployable adaptive SIGINT/IW sensor grid is 
different in applieation, but similar in approaeh to the wireless radar approaeh. Likewise, 
using beamforming teehniques for reeeption or eondueting offensive information 
operations of target signals of interest (SOIs) is different but similar in approaeh to 
beamforming for the purposes of data exfiltration. In this researeh, a novel approaeh for 
eondueting signals intelligenee from a distributed network of wireless nodes was 
developed. 

The primary objective of this research is enhancing RF collection in a specified 
target direction from a SIGINT/IW wireless sensor network. The proposed sensor 
network consists of a large number of small sensor nodes, K , distributed over an area, 

, to acquire information about signals of interest. These sensor nodes collaborate 
among themselves to form an ad-hoc network and disseminate the collected target 
information to a centralized controller node, e.g., UAV, UUV, etc. In this research, it is 
assumed that the clustering functions to create a sensor grid have been completed. It is 
assumed that perfect node location has been determined and that these data are available 
to the centralized controller node. 

To begin, this research explored the fundamental tenets of conducting 
SIGINT/IW. In order to enhance RF signal collection from this sensor network, these 
tenets were divided into broad categories. First is the determination of the frequency of 
interest of the target. Second is the synchronization of the individual nodes. Third is the 
determination of the direction of the target. Fourth, is the formation of the beam in the 
target direction. Two of these areas are specifically addressed in this research, categories 
three and four. The first category will be assumed to be accomplished by an 
implementation of the periodogram. The second category will be assumed to be 
accomplished through a “brute-force” method similar to the one presented by Jenn, et al. 
[21] for the wireless radar application. 
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To address categories three and four, two conflicting priorities were addressed. 
One is the time required to determine the target direction plus the time required to form 
the beams. The other is the energy consumption involved in developing these solutions. 
Two competing enhanced collection methodologies (ECM) for direction-finding within 
this wireless antenna array network were developed and analyzed. ECM-1 uses a 
combination of time difference of arrival (TDOA) and adaptive beamforming. ECM-2 
uses adaptive beamforming that performs a beamscan similar to phased-array radars. 

Analytical expressions for energy consumption and time to develop a solution, 
were derived to show that ECM-1 was far superior to that of ECM-2. ECM-1 was shown 
to be far superior to ECM-2 in both energy consumption and the time required to 
formulate a solution. 

In order to better focus on category four, two competing methods for forming the 
beams were developed. Method One selects a subset of sensors from a randomly 
distributed antenna array. This method uses data exclusively from these elements and 
increases the number of iterations until the desired beam is formed. Method Two creates 
a new subset of sensors, from a randomly distributed array, for each iteration, until the 
desired beam is formed. 

Analytical expressions were derived that proved that Method Two provides better 
power management across the sensor network. Additionally, Method Two was able to 
mitigate the grating lobes that are present in arrays where the antenna element spacing is 
greater than one half-wavelength. 

This dissertation is organized as follows. Chapter 11 will provide an overview of 
wireless sensor networks and lay the foundation for the SIGINT/IW wireless sensor 
system, to include characteristics and the proposed architecture. The chapter will also 
address some of the challenges in realizing this network. Chapter 111 will summarize 
previous work. The discussion will include digital array radar sensors, basic array 
antennas, communications architecture, wireless local oscillator element synchronization, 
and beamforming. Chapter IV will discuss an adaptive SIGINT/IW sensor network. The 
chapter will include a discussion on antenna theory and array properties, spectral 
estimation techniques and element synchronization. Chapter V will be a discussion on 
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subspace methods used for direetion of arrival determination. Chapter VI will diseuss the 
direetion finding teehniques of time differenee of arrival and beamforming. Chapter VII 
will diseuss the enhaneed eolleetion methodology developed in this researeh, and Chapter 
VIII will provide eonelusions and diseussions on future work. 
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II. OVERVIEW OF WIRELESS SENSOR NETWORKS 


As discussed in [2] and [68], wireless sensor networks eonsist of deviees that 
eombine the funetionality of sensing, eomputation and eommunieation into a single 
deviee eapable of self-organization and inter-deviee eonneetivity. Wireless sensor 
networks ean be used in a number of military and eivilian applieations. As seen in Figure 
1, these sensor nodes ean be deployed quiekly and left unattended by humans, whieh 
make them very attraetive for military applieations [2, 68]. Onee deployed, these 
SIGINT/IW sensor nodes have the ability to aequire information about SOIs and forward 
the information baek to a eentralized node for proeessing and relaying baek to a regional 
National Seeurity Agency/Central Seeurity Serviee (NSA/CSS) faeility, e.g., NSA/CSS 
Hawaii or a taetieal eolleetion platform. 

The ehapter is organized as follows. Seetion A will diseuss the proposed 
SIGINT/IW sensor network arehiteeture. Seetion B will provide an overview of 
SIGINT/IW sensor network oharaeteristies. Finally, Seetion C lays out some of the 
pereeived ehallenges in realizing this SIGINT/IW sensor network. 


Mobile Access Points 



Sensopi 


Figure I. Random deployment of sensor nodes by an unmanned aerial vehicle (UAV). 

[From Ref. 2] 
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A, SIGINT/IW SENSOR NETWORK ARCHITECTURE 

Several wireless sensor network arehitectures have been proposed in the 
literature. Figure 2 shows a general wireless sensor network arehitecture, which consists 
of a sensor field, a sink node and a satellite through which the user can access the sensor 
field [2], A two-tiered hierarchical clustering sensor network architecture is the one 
envisioned in this research. A central controller node (i.e., UAV, UUV, etc.) provides the 
link between the satellite and the sensor nodes, and each sensor node has the ability to 
route the data to the central controller node. The central controller node in turn sends the 
data to the user node via the satellite [73]. 

The central controller will be the reference antenna and is tasked with determining 
the frequency of the signal-of-interest, maintaining frequency, phase and data 
synchronization among the remaining nodes within the cluster. The network is comprised 
of K sensors spread over an area, . As seen in Figure 2 the target emitter lies far 
outside the sensor field. This is in contrast to most traditional sensor networks that 
monitor events that traverse through the sensor field. 



Figure 2. Proposed wireless SIGINT/IW sensor architecture. [After Ref 2] 

The requirements of sensor nodes dictate the need for an efficient system design. 
The sensors themselves will require much power, so minimizing the organizational 
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functions, and data exfiltration power requirements will be extremely important. This 
design eriterion also governs the intereonneetion of sensors. 


B, CHARACTERISTICS OF A SIGINT/IW SENSOR NETWORK 

As stated before it is advantageous to separate the proeess of eondueting 
SIGINT/IW into manageable seetions. The SIGINT/IW proeess was divided into four 
broad eategories. First is the determination of the frequeney of interest of the target. 
Seeond is the synehronization of the individual nodes. Third is the determination of the 
direetion of the target. Finally, is the formation of the beam in the target direetion. Two 
of these areas are speeifieally addressed in this researeh, eategories three and four. The 
first eategory will be assumed to be aeeomplished by an implementation of the 
Periodogram. The seeond eategory will be assumed to be aeeomplished through a “brute- 
foree” method similar to the one presented by , Jenn, et al. for the wireless radar 
applieation. 

The SIGINT/IW sensor network will be a eolleetion of deviees that onee deployed 
ean self-organize allowing for unattended operations. This self-organization will require 
the nodes to arrange themselves into a funetioning network. The network will be eentrally 
eontrolled. The nodes have the ability to eolleet radio frequeney information and 
eommunieate this information to the eentral eontroller. 

A SIGINT/IW sensor node will be assumed to possess the following 
distinguishing oharaeteristies. These distinguishing eharaeteristies will be the ability to 
demodulate and analog-to-digital eonvert intermediate frequeney (IF) signals, as well as 
the eapability to modulate and digital-to-analog eonvert IF for the purposes of eondueting 
offensive information operations. The sensors will have to have the ability to store large 
amounts of data required for the SOI, direetion finding, beamforming, ete [4, 57]. 

The storage eapaeity requirements will be mueh greater than those of typieal 
sensor networks in order to aeeount for ehannel aeeess and proeessing delay at eaeh 
sensor node, propagation delay between nodes, formulation of a given solution at the 
eentral proeessor and the ability of a third-party, i.e., an operator, to determine if 
sustained eolleetion is required [17]. This type of sensor network must operate as power 
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efficiently as practical to extend battery life. It will be assumed in this research that the 
sensors possess the capability of obtaining power from the surrounding environment, e.g., 
solar, sea, etc [6, 27, 65]. 

C. NETWORK CHALLENGES 

In order to realize the SIGINT/IW network, several network challenges will need 
to be addressed; these include, but are not limited to medium access control, localization, 
and synchronization techniques. Medium access control is the method through which the 
various nodes access the medium in order to relay data. Localization is the ability of a 
node to determine its physical location. Node synchronization is important in order to 
function as a coherent antenna array [6, 35]. 

1, Medium Access Control (MAC) 

The communication functionality of a sensor network node follows a layered 
protocol architecture. Figure 3 describes the typical layered protocol stack of a sensor 
network node. As laid out in [65], the application layer provides mechanisms for analog- 
to-digital conversion. The network layer is responsible for seamless transfer of 
information, and the data link layer provides fair access and is responsible for error-free 
transmission. The physical layer provides a means of sending and receiving a bit stream 
[64]. 



Figure 3. The protocol stack for a typical sensor network. [From Ref. 55] 
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Contention management is important in helping to avoid or recover from 
collisions, and minimize the overall penalty in terms of wasted network resources in 
shared-medium environments, such as contention-bus or wireless networks. Typically, 
these methods include random access, demand assignment, and fixed assignment [11], 

Due to the expected high sampling rates of the individual SIGINT/IW sensor 
nodes and the perishiability of the data collected, choosing or designing a MAC suitable 
for this type of network will be paramount. Each type of MAC scheme has advantages 
and disadvantages, but it is envisioned that a fixed assignment MAC will provide the best 
service. This idea is based on the notional fact that the data would be arriving at each 
node simultaneously and all of the data is important. All of the data is required by the 
central node in order to form the appropriate solution. 

2. Localization 

One of the primary objectives of deploying a SIGINT/IW sensor network is to 
detect the presence of SOIs and identify the positions of their emitters. This, in turn, 
poses a requirement on the sensor nodes to determine their own positions, thus some 
method for location discovery is needed. This section provides an overview of the 
relevant position location techniques and methods. More complete reviews and algorithm 
comparisons have been presented in [28, 38, 41, 44, 61]. Sensor networks require the 
extension of existing methods by finding ways to use measurements such as range or 
angle between pairs of unknown-location nodes. 

In a wireless environment, networked transceivers can obtain distance and angle 
measurements by one or a combination of the following ways: 


1. A sensor node’s location can be calculated directly and independently, if the sensor 
node is in range of multiple reference nodes. The sensors would calculate the 
estimated range by measuring the received signal strength (RSS) when a reference 
signal reaches the object [41]. The RSS can be measured as part of normal 
communications. Therefore, no additional energy or bandwidth is required. However, 
due to multipath reflections of the reference signal this technique can be 
unpredictable. 

2. The sensor node’s location can also be computed from the time-of-arrival. Time-of- 
arrival is simply the measure of the moment at which the signal first arrives at the 
receiver. The time delay between the sender and the receiver can be used to calculate 
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the angle of arrival or the distance when compared to a reference node This technique 
requires accurate measurement of the time-of-arrival of the line-of-sight (LOS) 
signal. This technique can also suffer from multipath and noise. 

3. Finally, a sensor node may use the angle of arrival to calculate position. The angle of 
arrival (AOA) technique requires accurate angle information, typically using phased 
antenna arrays with multiple antennas of known separation to perform angular 
calculation [28]. 

Due to the extreme requirement for perfect sensor location, timing 
synchronization and covertness, each node will be assumed to have some form of 
precision location capability such as a GPS receiver. This capability will ensure all nodes 
are able to be localized with a minimum of communications between sensors. This 
eliminates the requirement for beacon signals, thereby reducing the adversary’s ability to 
detect our network. It will be shown in Chapter III, Section G, subsection 2, how adaptive 
beamforming is affected by node location errors. 

3. Synchronization 

A persistent technological challenge arising from the development of the wireless 
opportunistic array concept has been the need to perform synchronization of the array 
elements to provide time and frequency references. Synchronization allows nodes the 
capability to determine their relative position when deployed randomly. To achieve data 
aggregation, the sensor must be able to precisely determine the instant in time at which an 
event occurred in order to recognize duplication [35, 65]. Previous research [14] has 
shown that control of the elements’ phase is possible via a wireless LO signal. In [28] 
two synchronization techniques were studied for a dynamic environment where the 
transmission paths will be changing and unpredictable. Generally, the techniques aim to 
get phase corrections for synchronization without explicit measurement of the elements’ 
position displacement. These two techniques will be discussed in more detail in Chapter 
III, Section E. 

This chapter provided an overview of the wireless SIGINT/IW sensor network. 
The architecture for this network was discussed to lay the foundation for the rest of the 
dissertation. Some of the general characteristics of the basic functions and missions of the 
wireless sensor network were covered. Finally, some of the networking challenges, not 
addressed in this research, but nevertheless are critical to networking, were examined to 
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give an appreciation for the complexity of this network. Now, a discussion of previous 
research relevant to the wireless sensor network will be offered. 
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III. SUMMARY OF PREVIOUS WORK 


This chapter summarizes some of the ongoing researeh relevant to the wireless 
SIGINT/IW sensor network. Seetions A - E are eomprised of wireless radar researeh 
eondueted by Jenn, et al. [20] The foeus in those research areas was on the sensor 
elements, networking and synehronization. 

Seetion F reeaps some of the time difference of arrival researeh eondueted by 
Loomis, et al. [29] Finally, Seetion G will diseuss the researeh into beamforming 
teehniques eondueted by Tummal, MeEaehen, Vineent, et al. [68] 


A, DIGITAL ARRAY RADAR SENSORS 

The opportunistie array researeh has led to the design of a digital transmit antenna 
with hard-wired array elements that operate from 2 to 2.5 GHz. The antenna was 
eonstrueted using the Analog Deviees AD8346EVAL quadrature modulator boards. 
Computer simulations using a genetie algorithm (GA), whieh has advantages for pattern 
formation when the array geometry is random or aperiodie, were used to verify formation 
of the radiation beam with randomly loeated elements. This part of the researeh 
demonstrated the viability of the transmit eomponent of the phased array [28]. 

In addition, the bandwidth eharaeteristies of the AD8346EVAL modulator board 
were investigated. Another eommereial produet, the Analog Deviees AD8347EVAL 
quadrature demodulator board was eonfigured to operate as a reeeive phase shifter. 
Referenee [28, 40] investigated the design of the eomplementary phased array reeeiver 
arehiteeture using the AD8347 demodulator. Different types of time-varying phase 
weights for a linear frequeney modulated signal were demonstrated for the transmit side 
to improve the phase distortion and inerease the operating bandwidth of the phased array 
[28]. 

While the previous researeh was foeused on the ereation of a wireless radar it is 
easy to see how the motivation is different but the technology is similar to the teehnology 
required for a SIGINT/IW sensor. For elarity the phrase SIGINT/IW has replaeed the 
phrase wireless radar to demonstrate how these nodes have potential applieations for use 
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in the wireless SIGINT/IW antenna network. A basie bloek diagram of the antenna array 
system arehiteeture is shown in Figure 4. This bloek diagram is an adaptation of the one 
presented in [28]. For elarity, only a single wireless module with a central controller is 
shown. The central controller should be capable of communicating wirelessly with 
hundreds or even thousands of array elements that are self-standing SIGINT/IW modules. 

For remote deployment applications, the central controller can be located in a 
UUV, UAV or a non-mobile central controller, while the array elements will be randomly 
distributed. 


Afiay ri«nft i 


'L 

to and Tmi ^ii:iL Rcvri^c 

Pluiw SyiiL 



I r.wsmMsjiW ^ 

/ z' 


/ / 
/ 

/ 


B\\nmforfnmg 
andloiitrolK'f Una 


LO t>ul 


‘i 


TdfgeleJ 

LOib 


WiiL'Ie'x'k PbLi 

to 

Llcrncm 


-11 


Y 


HcIciM Data 
Frcm 
l.lcmcnl 


f)lJltjl Itc.UMl'i’fllK't 

tnf Con(ri>i;« 


Figure 4. Network architecture for the SIGINT/IW. Only one sensor displayed. [After 

Ref 20] 

The general operational concept of the distributed array is as follows. At each 
array element, the received signal is downconverted to baseband after low-noise 
amplification, quantized (with the A/D) and the in-phase and quadrature data returned to 
the central controller for processing. Conversion directly to baseband will increase the 
sophistication of the sensor nodes, but it will drastically reduce the storage requirements 
at the nodes and it will reduce the network traffic, as well. The central controller 
computes the necessary control data - both phase and amplitude weights - for each 
element. In addition the central controller will maintain and deliver the appropriate signal 
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of interest (SOI) parameters. These are eombined with the loeal oseillator (LO) and 
synehronization signals and are passed wirelessly to all array elements. A detailed 
arehiteeture for eaeh SIGINT/IW sensor is shown in Figures 5 - 6. If offensive 
information operations are desired, the digital baseband signal is generated by the direet 
digital synthesizer (DBS), eonverted to analog (with the D/A), direetly up-eonverted to 
the operating band and power amplified. The LO referenee signal is distributed 
wirelessly. An aetive phasing teehnique is used to eompensate for element dynamie 
position ehanges and propagation ehannel ehanges. 


Wireless 

Data 
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Figure 5. System arehiteeture for the SIGINT/IW Tx/Rx module. [After Ref 20] 

The wireless opportunistie array eoneept demands very high data rate wireless 
eommunieation. Numerous self-standing SIGINT/IW sensors eontinuously eommunieate 
element loealization and synehronization signals, beam eontrol data, and digitized 
SIGINT signals wirelessly with the eentral signal proeessor [28]. A bloek diagram of the 
prototype SIGINT/IW sensor is shown in Figure 6 [28]. 
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Figure 6. Block diagram of the demonstration SIGINT/IW sensors. [After Ref. 21] 


B, BASIC ARRAY ANTENNA 

Jenn, et al. [28] created and tested a low-profde, broad-band U-slot microstrip 
patch antenna. The antenna was designed to operate in the upper VHF/lower UHF 
frequencies. A set of simple design procedures was proposed to provide approximate 
rules that result in a good “first-pass” design with prescribed characteristics that require 
minimal tuning. 


C. COMMUNICATION ARCHITECTURE 

In order to provide localization, pass synchronization, beamscanning and 
beamforming coefficients, other targeting data, etc., a wireless connection was 
established to integrate the elements throughout the opportunistic array. A proposed, 
communication architecture for the wirelessly networked opportunistic digital array radar 
(WNODAR) was presented in [28] and is shown in Figure 7. 
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Figure 7. Communication architecture for the opportunistic array. [From Ref. 28] 
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The data flow required to eonduet SIGINT/IW is analogous to if not the same as 
the data flow required to perform the funetions of the WNODAR. Periodieally eaeh 
sensor will send its position data to the eentral digital beamformer and eontroller via the 
position loeation data signal [20]. In the SIGINT/IW network this funetion is eontained 
within one system referred to as the eentral eontroller. When the eentral eontroller has 
determined the element loeations it begins to ealeulate the appropriate digital amplitude 
and phase weights for eaeh array element. The eentral eontroller then broadeasts this data 
via the waveform and eontrol data signal to all array elements [28]. The loeal oseillator 
information is required by the modulator and demodulator. Referenee eloek data is 
required for the DDS signals. This information has been eombined into a single 
waveform. A pulse train is transmitted from the eentralized eontroller and the pulse train 
envelope deteeted and used for timing [28]. The earrier ean be extraeted and used for the 
loeal oseillator. Eaeh Tx/Rx module will ineorporate hardware for performing the 
synehronization. The phase synehronization eontrol data is used to phase synehronize all 
the array elements [20]. The eentral eontroller will be used to establish a beamseanning 
pattern. When a signal is reeeived, the signal is demodulated and the SOI data is sent to 
the eentral eontroller for proeessing. [20] If transmitting then the amplitude and phase 
eorreeted waveform is then modulated, amplified and transmitted. A summary of the 
wireless data determined in [20] is shown in Table 1 [28]. 


Deseription 

From 

Data Rate 

Waveform eontrol 

Beam controller 

2 bit/s 

Synchronization control 

Beam controller 

1 bit/s for synchronization 
command 

1 bit/s for phase correction 
command 

Phase weights control 

Beam controller 

4 bit/s 

Received radar signals 

SIGINT/IW 

sensors 

16 bit X 2 xioo kS/s = 3,200,000 
bit/s = 3.2 Mb/s 


Table 1. Summary of wireless data requirements. [From Ref. 28] 


In [28], eontemporary wireless eommunieation systems were used for the 
prototype WNODAR sensors. A updated survey of the findings derived from [28] and 
[65] are presented in Table 2. It gives a brief overview of some eontemporary 
eommereial eonsumer wireless eommunieation systems. Current systems are not 
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sufficient to support both the high data rates and large numbers of nodes at the 
appropriate distanees required for the SIGINT/IW Network [28, 45, 65]. 


Table 


2 . 


System or 

Standard 

Frequeney 

Data Rate 

Standard 

Operational 

Range 

802.11a 

5 GHz 

54 Mbps 

100 m 

802.11g 

2.4 GHz 

100 m 

802.1 In 

2.4 GHz 

200 Mbps 

(typieal); 540 
Mbps (max) 

250 m (indoors) 

WUSB (UWB) 

3.1 to 10.6 GHz 

400 Mbps 

< 10 m 

Bluetooth 3.0 

6-9 GHz 

480 Mbps 

15 m 

802.15.4 

(Zigbee) 

2.4 GHz and 915 
MHz 

250 Kbps 

(max) 

70 - 300 m 


Performanee of eommereial eonsumer wireless communieation systems. 

Ref 28] 


After 


D, WIRELESS LOCAL OSCILLATOR (LO) SIGNAL DISTRIBUTION AND 

TRANSMISSION SYSTEM 

As previously stated the funetions of the WNODAR and the SIGINT/IW system 
are similar. One of the most fundamental and key eomponents of both systems will be the 
ability to maintain synehronization between the eentral eontroller and the sensors. One of 
the key tenets of this synehronization proeess will be the distribution of the loeal 
oseillator (LO). In the WNODAR projeet a method of wirelessly distributing the LO was 
explored. The teehnique from [21] is deseribed here to demonstrate a possible solution to 
the distribution of the LO that will be eneountered in the realization of the wireless 
SIGINT/IW antenna network. 

Distribution of the LO signal is required for eoherent array operations. As shown 
in [20], a wireless LO was sueeessfully demonstrated in a laboratory environment. A 
sinusoidal LO signal was transmitted, and then reeeived by two AD 8346 modulator 
boards operating at 2.4 GHz. During this laboratory experiment it was also shown that 
eontrol of the elements’ phase is possible via a wireless LO signal. 
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Several eritieal aspeets of the opportunistie array eoneept were demonstrated 
using COTS hardware in the 2.4 GHz band. Low eost, high-performanee modulators and 
demodulators are available at this frequeney. 

The LO referenee signal must be distributed to the modulators and demodulators. 
Synehronization of the LO referenee signal is neeessary to lower the sidelobe levels and 
improve the pointing aceuraey of the beam [28]. The LO referenee signal ean be 
distributed wirelessly; a synehronization eireuit is used to eontrol phase eorreetions due 
to element dynamie position ehanges and propagation channel changes. In the next 
section a more thorough review will be presented on the two synchronization methods 
presented in [28] to synchronize the elements in a dynamic environment. 


E. ELEMENT SYNCHRONIZATION 

As previously mentioned, synchronization of array elements to a common 
reference is required for the SIGINT/IW network to function as a coherent antenna array. 
Control of the elements’ phase is possible via a wireless LO signal, but in dynamic 
conditions the transmission paths will be changing and unpredictable. To overcome this, 
two different synchronization techniques were proposed in [28]. A brief overview of 
those techniques will be presented in this section. 

Synchronization of array elements in time and frequency ensures that the 
emissions from all elements converge coherently on the SOI, increasing average power 
and signal-to-noise ratio (SNR). In each SIGINT/IW sensor, the use of a DBS, 
demodulator and modulator will require precise phase-synchronization of multiple 
synthesized RF output signals to one another for coherent detection and integration. For 
our application, frequency is synchronized by a common wireless LO signal, and 
therefore the key focus is to provide time or phase synchronization. In [28], two 
techniques to perform element synchronization for the adaptive antenna array were 
investigated. 

The “brute force” synchronization technique [28] is a systematic adjustment of 
the array element phases. It is easy to implement with some hardware incorporated in 
each array element and in the central digital beamformer and controller. Figure 8 shows a 
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detailed diagram of the synehronization bloek that is required in eaeh of the SIGINT/IW 
sensors. Eaeh synehronization bloek eomprises a modem and eontroller eonneeted to a 
phase shifter and a switch. When the switch is positioned for synchronization operation 
(as shown) the LO signal is passed through a circulator, low-noise amplifier (ENA), 
phase shifter and then retransmitted back and compared to a reference signal at the 
central controller. Linder normal operation (switch opposite as shown), the EO signal is 
sent out to the modulator and demodulator for coherent beamforming. [28] 

At the start of the synchronizing cycle, the central controller sends out each 
element’s address in turn. When the element is selected, the switch is selected to the 
synchronization position, the received EO signal is shifted by (j)^ and sent back to the 
central controller. At the central controller, the EO signal from element n is compared 
with the received EO signal from the reference element. 


Coiwollet Data I 

Fcr Synchronization ~ 



Eigure 8. Diagram of a synchronization block for one element. [Erom Ref. 20] 


Eigure 9 shows the general concept to perform phase synchronization for element 
n . As described in [28], the EO signal from the central controller arrives at the synch 
block of each element at a different phase, , where is the wave number and r is 
the distance from the controller to the element n. One element is selected as the reference 

element and it receives the corresponding EO signal e . As stated in [28], the 
objective is to synchronize all the elements to the reference element by adjusting the 

phase shifter cj)^ to correct for the difference in path length -r^^- 
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Figure 9. Phase synehronization using the “brute force” technique. [From Ref. 20] 

Figure 10 shows the architecture for use with the beam tagging technique that was 
used in [28]. 



Figure 10. Phase synchronization using the “beam tagging” technique. [From Ref. 20] 

The performance of these techniques were tested for use with the SIGINT/IW 
network. These results will be presented in Chapter IV, Section D. 
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F. TIME DIFFERENCE OF ARRIVAL (TDOA) 

TDOA has been researehed extensively by Loomis, et al. [29] The primary foeus 
of this researeh has eentered around the determination of time delay between high-speed 
mobile reeeivers that are maintaining the same relative distanee to eaeh other. This 
researeh will leverage this previous researeh speeifieally in the areas of ealeulating the 
aetual time differenee utilizing the eross-ambiguity funetion and determining a line of 
bearing based on the estimated time differenee ealeulations using the Newton-Raphson 
Teehnique. This researeh will tailor these approaehes to the speeifie eharaeteristies of the 
SIGINT/IW network, speoifieally that the sensors are not moving. 


G. BEAMFORMING 

In [68], a three-tiered hierarehieal elustered sensor network arehiteeture was 
eonsidered. A simulation model was developed and implemented in MATLAB® to study 
the applieation of beamforming using sensor elusters to establish a eommunieation link 
with a UAV primarily utilizing a 2-D array. Different sensor node densities, position 
errors and sensor node failures were simulated to investigate their effeets on the antenna 
beam generated by the sensor eluster [68]. 

The results were eneouraging showing that the average sidelobe levels deereased 
while the maximum average power gain inereased when the number of sensor nodes 
inereased. In addition, the shape of the main lobe remained relatively unehanged when 
the number of sensor nodes ehanged.. However, it will be shown in this researeh that the 
mainlobe is highly frequeney dependent. Results also showed that the presenee of 
position errors and sensor node failures redueed the maximum average power gain of the 
antenna beam and inereased its sidelobe levels, but the shape of the main lobe remained 
relatively unehanged [68]. 

Follow-on studies [73] eonfirmed eonsisteney in the shape of the antenna main 
lobe under varying eonditions. Therefore the beam ean be steered towards the UAV as in 
the ease studied in [73] in an adaptive manner as the UAV flies over the sensor field. This 
implies that the energy of the antenna beam ean be direeted towards the UAV, resulting 
in higher signal-to-noise ratio, and, henee higher ehannel eapaeity. Additionally, the 
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maximum average power gain of the main lobe can be increased by increasing the 
number of nodes in the sensor cluster, thereby increasing the transmission range between 
the sensor cluster and the UAV [73]. 


Previously, a discussion was presented on how the node density affects the array. 
Now a discussion will be offered on how node element failures and node location errors 
impact the array. Geolocation depends on the knowledge of the elements’ position and 
this is even more crucial in digital beamforming. The individual elements are placed in 
open, available areas and the positions will be random. This fact must be taken into 
account in the signal processing to avoid degradation in the sidelobes, gain and beam 
pointing [28]. 


Now a presentation will be given, which analytically describes the effects of 
element failures and position errors. 

1. Sensor Failures 

Element failures will be a reality in sensor networks that will need to be 
mitigated. An analytic explanation articulated in [27, 68] will now be used. The fractional 
loss in the main lobe gain for a two-dimensional MxN uniformly displaced array 


with equal weights due to element failures is given by [27, 68] 

^MW " 


h = 


MN 


( 1 ) 


where M' N' is the number of operating elements and MN is the total number of 
elements. Since element failures are random, the operating probability of any element 
is given by 

^MN'^ 


Po= lim 

MN^oo 


( 2 ) 


7 


■V MN 

Hence, the expected value of the main lobe gain, relative to an error-free array, is 
reduced by a factor of [28]. Figure 11 shows the comparison between the 


performance a two-dimensional array with the elements separated by an equal distance 
and a two-dimensional array with elements uniformly distributed within a given area. The 
simulation used a sub-array of MxN = 250 nodes uniformly distributed within an area, 

= 25m^. 
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Figure 11. Fractional Loss in the gain of the mainbeam resulting from node element 

failures. 

The simulated results were slightly poorer than that of the analytical prediction. 
However, the results did consistently follow the curve as seen in Figure 11. The results 
demonstrate how robust the beamforming algorithm of [73] is in the presence of node 
failures. The antenna array was able to sustain the loss of up to 40% of its elements 
before significant reduction in the gain of the mainbeam. What is interesting to note is 
that even though the mainbeam in all cases was reduced; the sidelobes were reduced as 
well. This allows the antenna array to maintain a direction finding capability even in 
these adverse conditions. In the next section a discussion on the effects of node location 
errors will be presented. 

2, Sensor Location Uncertainty 

The fractional loss in the main lobe gain for a two-dimensional MxN array 
with a equal distance between elements and equal weights is given by [27, 73] 
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where is the varianee of the phase error, assuming gaussian phase errors [73]. As 

per [57], the varianee of the phase error is related to the variance of the node location 
error through the relationship shown in Equation (4). 



f 

V 



y 


(4) 


It can be seen in Figure 12 how the array with uniformly distributed elements 
within a given area using adaptive beamforming out performs the two-dimensional array 
that has equal spacing between elements. 


Fraction Loss of Main Lobe Gain (dB) 



Figure 12. Fractional Foss resulting from node location errors. 


Just like in the case of node element failures the fractional loss is not dependent 
on frequency. Figures 11 and 12 show the robustness of the FMS adaptive beamforming 
algorithm of [73]. This robustness is the exact properties desirable for applications within 
the SIGINT/IW network. 
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The research conducted in [68] was primarily focused on the exfiltration of data. 
The focus of adaptive beamforming in this research will be for the purposes of enhancing 
collection and offensive information operations utilizing a randomly distributed sensor 
array. Now a more thorough discussion will be presented on adaptive networks and 
direction finding. 

This chapter summarized some of the previous work that has been conducted in 
areas germane to this dissertation. The chapter offered insight into research being 
conducted into the individual RF signal collection sensors, both the receivers and the 
associated antenna and an associated synchronization methodology. This chapter also 
provided discussions of research being conducted in the area of a proposed 
communications architecture. Finally, the chapter looked at time difference of arrival and 
Beamforming. This dissertation will use the research conducted in both of these areas to 
develop the proposed SIGINT/IW wireless antenna system and the enhanced collection 
methodology. 
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IV. ADAPTIVE SIGINT/IW SENSOR NETWORK 


This chapter will provide an overview into some of the basies of eondueting RF 
signals eolleetion. First, a diseussion on basie antenna theory and antenna array properties 
will be provided. The eoneept of speetral estimation will be diseussed in order to provide 
the reader with a general understanding into how the proposed SIGINT/IW system will 
determine the frequeney being transmtted from the target. Finally, a diseussion will be 
offered on the synehronization of the individual sensor nodes. 


A, ARRAY THEORY 

It will be assumed that a single element will have an isotropie radiation pattern, 
whieh means the radiation is distributed equally in all direetions. Therefore, it has low 
values of direetivity and gain [3, 14]. Geoloeation systems require very direetive 
eharaeteristies and high gain to meet the demand of high resolution. Two methods ean be 
used to realize these requirements. As diseussed in referenees [3, 14], the first method is 
to inerease the eleetrieal size of the antenna. The seeond is to form an array by 
assembling the elements in a speeified pattern, but keeping the original size of the 
individual elements [3, 14]. 

Sinee the SIGINT/IW nodes are envisioned to be low power and minimally 
eapable, the network needs to take advantage of the numerous elements present within 
the array. An antenna array ean take advantage of the multi-elements and known 
separations between elements to eombine the outputs and affeet the antenna array pattern. 
This pattern will quite often differ from the pattern of the individual elements. An array 
has the ability to aehieve the direetional performanee of larger antennas. By varying the 
phase and amplitude of the individual element outputs before eombining, the overall 
array pattern ean be steered in the desired user's direetion without physieally moving any 
of the individual elements [3, 14]. 
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By the principle of pattern multiplication, the overall radiation pattern 
is found as the product of the individual element radiation patterns ^[<x),6,(/>) with the 
array factor 'F ( 

The array factor is determined by the positions of the elements relative to the 
central controller, as well as the phase and amplitude levels of the received signal. The 
array factor, and thus the overall pattern of the array, can be continuously scanned or 
adapted by adjusting the relative phase and amplitude levels between the elements [7, 14, 
60], 

A linear array is one in which the elements are aligned along a straight line. If all 
of the elements lie in a plane, the array is a planar array [7, 14, 60]. Examples include the 
linear array, a circular array, and arbitrarily shaped planar arrays. Due to its simplicity, 
the uniformly spaced linear array shown in Figure 13 will be used to show the phase 
relationship between the received signal at the different elements for a given angle of 
arrival (j)^ measured along the main axis of the array. Since the signal present at element 

two has traveled a distance d cos(^/i^) farther than the signal present at element one, its 

phase will lag behind that of element one by (id cos(^/i'J where jd = In ! X is the phase 

propagation factor or wave number. This relation holds for a narrowband signal. A 
narrow band signal is defined as a signal whose modulating signal bandwidth is much 
less than the modulated signal [7, 14, 60]. 

The narrowband assumption allows us to assume that the only difference between 
the signal present at element two and element one is the phase shift induced by the extra 
distance traveled and is not significantly affected by the modulation during this time. An 
example in [60] demonstrates this concept. The complex envelope of a narrowband signal 
is shown in Equation (6) [7, 14, 60] 
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= ( 6 ) 
where A(^tJ is the signal amplitude, /(t) is an arbitrary phase shift and cOg is the center 
frequency. Assuming isotropic elements, <^(cOg) = l, and taking the received signal at 
element one as the reference, the received signals h, (t) for a uniform linear array with 
element spacing d can be represented in matrix form as 




1 

^-7Afl'cos(4) 

^-y/32rfeos(^,) 




(7) 


where N is the number of elements, (j)^ is the angle of arrival, and is the array 


factor. For simplicity, the frequency dependence of 'F has been dropped for the 

narrowband case, and the elevation angle 0 is assumed to be zero relative to the 
element's Foresight. The array factor must be carefully measured to calibrate the array for 
direction finding experiments [7]. 
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Now this concept will be extended to an Mxl array with isotropie elements. Three 
areas in partieular will be discussed. Two areas involve the baseline spacing, d , grating 
lobes and mutual coupling. The third area involves eenter frequency jitter, which can 
cause beam squinting. 

1, Grating Lobes 

The spaeing between antenna elements is critical in the design of a phased array. 
This is beeause when the spacing exceeds certain eritical values, grating lobes will occur 
[31]. Grating lobes are additional sidelobes that have the same amplitude as the main 
beam. It is normally undesirable to have grating lobes, and in order to avoid them, the 
element spaeing is required to meet the condition given by 


where \ is the wavelength of the highest operating frequency and (j)^ is the sean angle. 

Based on this eondition, when scanning to endfire at (j)^ = ±90° the first grating 
lobe will occur when the spacing is greater than \ I2 . For a comparison see Figures 14, 
15 and 16. Figure 14 has \ spacing. Figure 15 has 3>\ / 2 and Figure 16 has 
[14]. 

It can be observed in Figures 14, 15 and 16, for the uniform linear spaced array 
used in the above caleulations, that the antenna pattern is changed drastically as the 
spacing for the antenna is ehanged. This direetly affects the amount of signal and noise 
that is eollected and further analyzed. As the number of main lobes increases, the signal 
and noise, are no longer uncorrelated. As will be shown later too much correlation exists 
between the signal and noise for the classical subspace algorithms like MUSIC and 
Maximum Likelihood Method to separate them into individual components. Other signal 
processing methods will need to be employed to mitigate these grating lobes. 
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Figure 14. Array pattern for a 10 element Uniform Linear Array with half-wavelength 

spacing of the elements at a frequency of 157 MHz. 
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Figure 15. Antenna pattern for a 10 Element ULA with one and half-wavelength spacing 

at a frequency of 157 MHz. 


31 



I tmwni 2A «l 1S7 Wi/ 



•■ao 


12 * 


•tM 


136 


•tM 


60 

lao ^ 


Figure 16. Antenna pattern for a 10 Element ULA with two times the wavelength spacing 

at a frequency of 157 MHz. 

The element spacing used in the above simulations was not a realistic spacing 
given the random nature of deployment of our sensors. The sensors in the SIGINT/IW 
network on average would most likely be much further apart. For a more realistic spacing 
of one sensor every 10 meters one can see in Figure 17 that the situation is only getting 
worse. 
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Figure 17. Array pattern of 10 element Uniform Finear Array with 5 times the 
wavelength spacing of the elements at a frequency of 157 MHz. 
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A two dimensional planar array is a better representation for the SIGINT/IW 
sensor network. As seen in Figures 18 - 20, the situation is only marginally better. The 
simulations indicate that the beamwidth gets narrower, but the number of mainlobes 
remains consistent. This is consistent with the fact that as the number of elements 
increase so does the gain and as the gain increases the beamwidth gets narrower assuming 
that the effective aperture continues to increase as the number of nodes increases. 
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Figure 18. 
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Antenna pattern for a 20 Element 2D Uniform Array with half-wavelength 
spacing at a frequency of 157 MHz. 
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Figure 19. Antenna pattern for a 20 Element 2D Uniform Array with one and half¬ 
wavelength spacing at a frequency of 157 MHz. 
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Figure 20. Antenna pattern for a 20 Element 2D Uniform Array two times the 

wavelength spacing at a frequency of 157 MHz. 
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As stated before the elements will on average not be spaced that closely together. 
As seen in Figure 21a more realistic spacing of about 10 meters clearly shows that the 
number of mainlobes continues to increase. 


20 undorm 20 ivrivy wifi 2X tpAnng 



Figure 21. Antenna pattern for a 20 Element 2D Uniform Array with five times the 

wavelength spacing at a frequency of 157 MHz. 

Increasing the number of elements decreases the width of these individual lobes, 
but does not decrease them in number as can be seen in Figure 22. 
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Figure 22. Antenna pattern for a 40 Element 2D Uniform Array with five times the 

wavelength spacing at a frequency of 157 MHz. 


Research suggests that this problem only gets worse as the frequency increases, as 
can be seen in Figures 23 and 24. 
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Figure 23. Antenna pattern for a 20 Element 2D Uniform Array with 10 meter element 

spacing at a frequency of 800 MHz. 
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Figure 24. Antenna pattern for a 20 Element 2D Uniform Array with 10 meter element 

spacing at a frequency of 2.4 GHz. 

2, Mutual Coupling 

Mutual coupling occurs because the antenna elements within the array interact 
with each other. This interaction between elements results in an impedance change as 
seen by each element, which in turn affects the current magnitude, phase and distribution 
on other neighboring elements [53, 60]. 

In general, the element spacing should be designed to avoid grating lobes and 
reduce the adverse effects of mutual coupling. According to [60], the spacing is 
recommended to be between 2-/3 and XU. 

Due to the random distribution of our sensors, the distances between antenna 
elements will be such that mutual coupling should not be a problem. However, grating 
lobes will be something that will have to be mitigated as node densities get less. The 
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sensor grid will in some cases be spaced at distances larger than XI2. The choice of 
spectral estimation and direction-of-arrival technique will need to be capable of coping 
with this reality. 

3, Beam Squint 

Phase shifters are designed to adjust phase shifts that may occur during 
transmission to a desired center frequency and, if the signal that is being received is 

not at /q as may be the case for wideband signals, an effect called “beam squinting” will 

occur. An example taken from [40], of beam squinting is shown in Figure 25, where the 
beam is first pointed to 20 degrees and the frequency is then changed to 1.7 GHz, .8 GHz 
and 2.5 GHz. It can be seen that the scan angle of the main beam decreases for 
frequencies higher than the center frequency and increases for frequencies lower than the 
center frequency. The beamwidth also becomes narrower at higher frequencies and wider 
for lower frequencies. It can be seen that when the operating frequency changes, then the 
main beam will steer off from the desired direction, changing both in width and direction 
with frequency [49]. As stated in [40, 49], the beam squint is greater for a larger array 
aperture than a smaller array aperture. 



Figure 25. Beam patterns for a phased array at 0.8 GHz, 1.7 GHz and 2.5 GHz when 
phase shifters are set to steer beam to 20 degrees at 1.7 GHz. [From Ref 40] 
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Figure 26. Array pattern for a phased array at 157 MHz, 800 MHz and 2.4 GHz using 

adaptive beamforming. 

o 

In Figure 26, where the beam is first pointed to 45 and the frequency is then 
changed to 157 MHz, 800 MHz and 2.4 GHz. It can be noted that by using an adaptive 
beamforming algorithm that the scan angle does not drift, however the beamwidths are 
still different. The beamwidth also becomes narrower at higher frequencies and wider for 
lower frequencies. 

The above examples were for the narrowband case, but references [40, 49] 
suggest that it is possible for a phased array to achieve wideband performance by 
changing the settings of the phase shifters whenever the frequency of a signal with 
narrow instantaneous bandwidth is changed. This is equivalent to radiating multiple 
narrowband signals one at a time over a wide range of frequencies by adjusting the 
settings of the phase shifters. 
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B, ARRAY PROPERTIES 

1. Linear Array 

Now a more in depth discussion of array properties will be presented. Continuing 
with the narrowband example, the array depicted in Figure 13 has been extended to a 
Mxl array shown in Figure 27. The Mxl array shown in Figure 27 is located in the far 
field of a point source. The point source (target signal) arrives at an elevation angle with 
respect to the array normal (i.e., z-axis). For convenience, the reference (central 
controller) is taken at the origin, and it is assumed that the wavefront arrives at the 
element seconds before it reaches the origin. Hence, the signal arriving at the element 
leads the signal arriving at the origin by [73] 




c 


(9) 


where c is the speed of light, 0^ is the direction of arrival and is the x-coordinate of 
the element. 
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Figure 27. AnMxl omni-directional antenna array. [After Refs. 60,71] 


For wideband signals, the parameters in Equation (9) vary with frequency, so now 
the array factor is given by 

M 

^>{ 0 ,( 0 ) = , ( 10 ) 

m=\ 

where (X) = 2nf is the center frequency of the signal and is the element weight 
amplitude [49]. 

If the frequency of operation changes, then a new set of phase shift settings is 
required to keep the beam location fixed. This means the first exponent in Equation (10) 
must change with frequency in the same way that the second exponent changes with 
frequency. 
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It can be seen from Equation (9) that ) is determined by the angle of arrival 

of the desired signal indieated by 6^ and the element’s position indieated by . For any 

arbitrary angle 6, the signal arriving at the m‘^ element leads the signal arriving at the 
origin by 


>.( 6 ) 


COS 6 

c 


( 11 ) 


where, continuing with the narrowband example, the co dependence has been dropped. 
The mutual eoupling effects have been ignored. The array factor of the is obtained by 
adding all the array element outputs together giving [60,73] 

M 

vp i^e) = ^ ^ J2) 

m=l 

where and magnitude and phase of the current indueed on the 

element. Together, and form what is eommonly known as the eomplex 

weights. Notice that the phase reversal in the term jg required in order to ereate 

a maximum value of 'E (^) at 0 = 6^. Henee, the maximum value of of the magnitude 

'E (^) oeeurs at 0 = 6^ and the main lobe points towards 6^. Equation (8) can be written 
in the form given by [60, 73] 

M 

m=l 

where jg complex weight applied to the m'* element. The far-field 

power gain is given by 


\F{0t = 

where 7^(6*) is the array pattern with isotropie elements defined earlier. If the Mx\ array 
has the main lobe is pointing at 6 = 6^, the maximum array factor is computed to be 


M 

z 




-jP{^ 






(14) 
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(15) 


MM M 


'r('?„) = Z»’-,£'^'’‘“-‘“’ = Z4 


1,1. =M, 


using Equation (13) while the power gain is ealeulated using Equation (14). Therefore, 
the power gain of a uniform array ean be inereased by inereasing the number of elements. 
This implies that more sensor nodes can be used within a sensor cluster to produce an 
antenna beam with higher power gain and hence increase the detection range and 
offensive information operations capability between the sensor network and the target. 

To illustrate the main characteristics of the antenna beam, a plot taken directly 
from [73] is utilized. The plot of the normalized power gain (in dB) versus 0 (in 

degrees) is shown in Eigure 28 and was generated by a IQx 1 uniformly excited array. The 
array consists of isotropic elements with fixed element-element spacing d = XU, and the 
main lobe, as shown in Eigure 28 (a) points in a direction perpendicular to the array axis 
(i.e., x-axis) and can be categorized as a broadside radiation pattern. As shown in Eig. 54 
(b), the 3-dB (or half power) beamwidth is about 10.2°, and the highest sidelobe 

level 5*^^^ is about 13.2 dB below the main lobe level. 

The 3-dB beamwidth of a uniformly excited broadside array is given by [60, 73] 

= 0 . 866 ^, (16) 
Md 

where a long array (i.e., Md n X) has been assumed. The 3-dB beamwidth was 
calculated to be 9.92° using Equation (16), showing reasonable agreement with the 10.2° 
obtained from the plot. It can be seen from Equation (16) that 6 * 3 ^^ can be controlled by 

adjusting M . Eor example, when M is increased by a factor of 4", reduced by 

the same factor, resulting in a narrower beam and vice versa. This allows 6 * 3 ^^ to be 
adjusted by changing the number of sensor nodes required for beamforming during the 
detection and tracking of UAV. Eor example, a large is usually required to create a 

large detection region while a small 6^^^ is preferred to localize & track the target and 

perform sustained collection. What is assumed is that when the number of elements 

increase so does the effective aperture. It will be shown in Chapter VI, Section B, 
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Subsection 3 that when the physical area is held constant and the number of nodes 
(elements) is increased the beamwidth actually gets broader and not narrower implying 
less gain. This is further discussed in sub-section 4. 




Figure 28. Normalized power gain of the beam generated by a IQx 1 uniformly excited 
square array with isotropic elements, fixed element spacing oid = X/2, and 9a = 0° 
(a) Polar plot and (b) Rectangular plot in the X-Z plane indicating the 3-dB 
beamwidth of 10.2° and the highest sidelobe level of -13.2 dB. [From Ref. 73] 


The array factor and the power gain of a one-dimensional array can be extended 
to a two-dimensional array as presented in the next section. 

2, Planar Array 

Figure 29 shows an MxN omni-directional antenna array located in the far field 
of a point source. The received signal arrives at an elevation angle, 0 with respect to the 
z-axis and an azimuth angle, (j) of with respect to the x-axis. The central controller is 
taken at the origin, and it is assumed that the received signal arrives at the element 
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seconds before it reaches the central controller. Therefore, the signal arriving at the 
{^m,nY element leads the signal arriving at the origin by [60, 73] 

. (n .. ^_^«,«sin6»,cos^^io+y„„sin6»^sin^o 

’ (i/) 

where c is the speed of light and is the x-coordinate and y-coordinate of the 
element. 



Figure 29. AnM^N omni-directional antenna array. [After Ref 73] 


For any set of arbitrary angles, , the signal arriving at the element 

leads the signal arriving at the origin by 

sin 6 cos (h + sin 0 sin S 

mn T J mn T 




The two-dimensional array factor is given by 

M N 






-jP{x^^ sin^) sin^cos^+T^„ sin^sin^) 


(18) 


, (19) 


=1 n=\ 


without considering the effects of mutual coupling among the array elements. [60] 
Equation (19) can be rewritten as 


M N 




• JtU- sin^cos^+T„„ sin^sin^ 


( 20 ) 


m=l n=\ 
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where (5 is the wavenumber given by In! X and 

the eomplex weight applied to the element. The maximum value of F{0,4>) 

oeeurs at {^6,(/)) = {^9^,(/)q), and the main lobe points towards {0^,(/)q), The radiation 
pattern with isotropie array elements is defined as 






M N 

II 

ff 2 =l n=\ 


* sin^cos^+Tmw sin^sin^ 


w e 

mn 


( 21 ) 


Consider an MxN uniformly exeited array with isotropie elements and the main 
lobe pointing at 0^=(/>q= 0° The maximum array faetor and the power gain are () 

and (^NMy, respectively. This result is similar to the linear array. The same conclusions 

can be drawn that the power gain of a planar array can be increased by increasing the 
number of elements with the assumption that the effective aperture increases. 

The following example taken from [68] illustrates the affects of different array 
sizes on the beamwidth of the beam. The beam pattern formed by 3x3, 5x5 and 7x7 
arrays are shown in Figure 56 in which plots were generated for (j) = 45°. These 
uniformly excited square arrays consist of isotropic elements with a fixed element- 
element spacing d = XU. The main lobe for all the three arrays is assumed to point at 

K.A) = (o'.«-) as illustrated in Figure 30(a). Figure 30(b) shows that 6],^^ increases 

from 15° to 37° when the number of array elements decreases from 7x7 to 3x3. Therefore, 
the antenna beam of a uniformly excited array with isotropic elements and fixed element- 
element spacing becomes broader when the number of elements decreases. This 
observation is consistent with the conclusion in the previous section for the one¬ 
dimensional case. Section 4 will investigate the beamwidth affects from a randomly 
distributed network. 

In summary, the array factor and power gain equations for the planar array were 
presented in this section, and it was concluded that they were affected by the azimuth of 
the target signal, the node density, the node spacing, and the complex weights Wmn- The 
main characteristics of the antenna beam, such as 3-dB beamwidth and highest sidelobe 
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level, were introduced. The effects of the number of array elements on the antenna power 
gain and the 3-dB beamwidth were also discussed. 






7x1 •fray 5x5array ^x^arraiy 2Dj>lo»4gcncrat«l 

Figure 30. Normalized power gain of the beams generated by 7x7, 5x5 and 3x3 

uniformly excited square arrays with isotropic elements, fixed element spacing of 

d =XI2, and (6 *^,) = (o°,45°) : (a) Polar plots and (b) X-Y plots showing the 
decrease in 3-dB beam-width as the number of elements increase.[ From Ref 73] 

3, Random Array 

The SIGINT/IW network will consist of SIGINT/IW sensor nodes that are 
deployed randomly over an area of interest, resulting in random positioning of the sensor 
nodes. Subsequently, these sensor nodes form clusters and their locations are determined 
using location discovery techniques and are reported back to the central controller. The 
random position of the sensor nodes within the sensor cluster gives rise to what is 


48 









































referred to as a random array. In addition, the SIGINT/IW antenna elements will in some 
eases not be loeated within a half-wavelength from each other. Therefore an algorithm 
capable of mitigating the grating lobes discussed in section A, subsection 1 will have to 
be employed. 

The SIGINT/IW sensor network will be a random array antenna that will 
adaptively control the pattern in order to enhance collection and/or facilitate offensive 
information operations in the direction of the target. Once the target has been localized 
the central controller can enhance the gain toward the desired signals and can null toward 
undesired signals. Correct control of beams and nulls is indispensable for the adaptive 
antenna. The pattern of the array will be controlled by dynamically varying the phase and 
amplitude of the received signal of each element. There are a variety of weight control 
algorithms for adaptive antennas. This research will build on the LMS algorithm research 
conducted in [73]. 

Continuing with the narrowband example from Section 1, each element output 
s. (t) is multiplied by a complex weight w.*, modifying the phase and amplitude relation 

between the branches, and summed to give h(t). 


'P (6>) = [wi* w/ W3* 


1 



^-jPdco%(e) 

^-jP2dcos(9) 


b{t) = w^s{t), (22) 


y^-jP(N-\)d^os(e) ^ 

Since the array factor for the overall antenna pattern is dependent on the phase 
and amplitude relationship between the branches and the weight vector modifies the 
phase and amplitude relationship, the overall array pattern can be continually modified by 
adjusting the weight vector [14, 60]. The overall pattern for the array is given by 
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(23) 


F{(o,6,(j)) = 6',^z5)| , 

by changing the phase in each element the network can steer the beam virtually 
instantaneously. 

4, Controlling the Size of the Beamwidth 

Understanding how to control the beamwidth of the beams will be eritical. In a 
seareh mode a wider beam, i.e. larger detection region is preferable. When sustained 
eolleetion and/or offensive operations is desirable then a more directive beam with higher 
gain in the target direction will be desirable. The equations for a two-dimensional planar 
array taken from [3] will be used as a guide for comparison with the randomly configured 
wireless antenna array. 

The beamwidth for a planar array is determined by using Equations (24) and (25). 
In order to determine the respective beamwidth the array is separated into x and y linear 
arrays. These equations for the respective scan angles for both elevation and azimuth are 
shown in Equation (26) [3]. 


0 


elevation 


cos" (6»^) \e^-^ cos" sin" (^)_ 


(24) 


0 


azimuth 


1 


sin" cos" )_ 


(25) 


where is the desired azimuth scan angle and 6^ is the desired elevation angle. It can 

be seen how the beamwidth in the azimuth and elevation direetions impaets the gain of 
the beam. [3] 


6^ =eos ' 
= cos ' 


cos 6*^ -0.443 
cos 6*^ +0.443 


{J^Length + d^ j 


yi^Length + d^ J 


( 26 ) 
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where d = Length / N , Length is the length of the uniform linear array and N is the 
number of nodes along the length of the uniform linear array. Now shown below for the 
azimuth [3]. 


=cos 

/ Xq 

^/i^^=eos 

where d = Length / N, Length is the length of the uniform linear array and N is the 
number of nodes along the length of the uniform linear array. It will be shown that as the 
number of elements inereases for a specific area, i.e. the node density increases, the 
beamwidth reaches an upper-bound for a given wavelength [3]. 

Likewise it can be seen from Equations (26) and (27) that as the length of each x 
andy increases for a given number of nodes that the beamwidths get narrower. Equation 

(28) is for the solid beam of the array using the azimuth and elevation calculations from 
Equations (24) and (25). Equation (29) is for the directivity, gain in a particular direction, 
which is based on the solid angle of the array. It can be seen from Equations (28) and 

(29) that as the beamwidths get narrower the directivity goes up increasing the gain in 
that particular direction and vice-versa when the beamwidths increase the directivity in a 
particular direction goes down. 

Armed with those two vital facts the central controller can now control the width 
of the beam by selecting areas of the overall grid that contains a sufficient node density in 
order to produce the appropriate beam size. This feat can be accomplished simply by 
choosing nodes that are contained within a specific area. 


cos^zig -0.443 

L 

' ^cos^zig +0.443 


^ zl ^ 


(L + t/) 


V 

^ zl ^ 


(L + (7) 


(27) 


^ A ^ elevation ^ azimuth 


(28) 


DU 
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(29) 



The inverse relationship between directivity and beamwidth can be observed by 
analyzing Equation (29). Figures 31 (a) - (c) show this inverse relationship between 
directivity and beamwidth for three frequencies: 157 MHz, 800 MHz, and 2.4 GHz. The 
three frequencies were distributed over different areas, A 2 = 25m ^, = 9m ^, A^ = 4m ^, 

respectively to account for the differences in wavelength. The simulations ignore any 
mutual coupling effects that would occur as the baseline, d , became less than /I / 3. This 
relationship is important because for a search mode a larger detection region, i.e. larger 
beamwidth, is desirable. The ability of the central controller to control the beamwidth 
will be an important operational constraint. 



(a) 157 MHz 
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(b) 800 MHz. 



(c) 2.4 GHz 

Figure 31. The inverse relationship between direetivity and beamwidth. 250 nodes were 
distributed over three different areas depending on the frequency tested, (a) Area, 
A 2 = 25 for frequency =157 MHz; (b) = 9 for frequency = 800 MHz; 

and (c) = 4 m , for frequency = 2.4 GHz. 
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For a given SIGINT/IW network with K sensors spread over an area, , the 
average distance between the sensor nodes will increase as the central controller creates 
an array containing A sensors contained with an increasing area, , to begin 

beamforming. As can be seen in equations (24) - (27) as the average distance increases 
between sensor nodes this causes a decrease in node density and a resulting decrease in 
the beamwidth. This relationship is illustrated in Figures 32-34. 

B»3mwidth w dislanc» betw»«n nod»s. 157 MHz 



Figure 32. The relationship between beamwidth and node density. The Average distance 
between nodes increases as node density decreases. Nodes are contained with an 

Area, A 2 = 25 square meter area. 
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Beamwidih v? Number of nodes, 157 MHz 



distance between sensors ■ (m) 


Figure 33. The relationship between beamwidth and node density. The Average distance 
between nodes increases as node density decreases. Nodes contained within an 
Area, Jj = 2500 square meter area. 


Beamwidth vs distance- between nodes, 157 MHz 



Figure 34. The relationship between beamwidth and node density. The Average distance 
between nodes increases as node density decreases. Nodes contained within an 
Area, J 2 = 250000 square meter area. 
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Simulations in Chapter VI, Section B will be used to illustrate how the beamwidth 
of the mainbeam is affected by the central controller. The central controller can affect the 
mainbeam by changing the node density through its selection process. Specifically, as the 
number of nodes is held constant and the central controller selects from an array of nodes, 
A, that are dispersed within a smaller area, , the beamwidth of the mainbeam can be 

widened for search. In effect the central controller is simply selecting the nodes that are 
closer in distance to itself. Likewise, as the central controller selects from an array of 
nodes, A, that are dispersed within a larger area, , the beamwidth of the mainbeam 
can be narrowed for sustained collection. In effect the central controller is selecting nodes 
that are located further from its location. 

In the above case the effective area remains constant and SIGINT/IW sensors 
were added in effect increasing the node density. Since the node density increases and the 
effective area remains unchanged the baseline distance, d , between nodes decreases. The 
larger the baseline, d , between nodes results in a narrower beamwidth which is similar to 
the previous sections where the nodes become spread over a greater area, A ^, then the 
beamwidth gets narrower. Therefore as the node density increases the baseline, d , 
becomes smaller resulting in a broader beam, which can be seen in Figures 32-34. This is 
illustrated quite well in Equations (26) - (27) that as the number of nodes is increased. 
Equations (25) - (26) asymptotically reach a constant value. This will be illustrated in 
Chapter VI, Section B. 

It can be assumed that a reasonable node density will be required in order to do 
adaptive beamforming. As will be shown in Chapter IV, Section B this node density is 
highly dependent on the frequency of the signal-of-interest. This highlights the 
limitations of this type of network in using beamforming to enhance collection at higher 
frequencies. The node densities required for enhanced collection would be high, which in 
turn would require the sensors to be small to maintain some level of covertness. This 
relationship will be further discussed in Chapter IV, Section B. 

It will be shown in Chapter VI that the beamwidth of the beam can be controlled 

by changing the subset selection area of the wireless SIGINT/IW antenna network. Eor 

our sensor grid, the central controller selects the number of array elements needed to 
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form the desired beam. The central controller in effect creates an array, which is merely a 
subset of the total array. For example, consider a sensor grid that contains K elements 
randomly spread over a total area, A-^. The central controller will select nodes from an 

area containing A elements randomly spread throughout the area, where A<K and 
A 2 < A^. Of the nodes within area, A 2 , F nodes will be chosen, where PU A. The P 
node subset will form the basis for the array factor. The area A 2 , in effect becomes the 
effective aperture of the antenna grid. Based on the size of the area, A 2 , the beamwidth of 

the beam can be tailored to specific applications, specifically, a wide beam for detection 
and a narrow beam for localization and/or collection. 

C. SPECTRAL ESTIMATION 

The periodogram and the Multiple Signal Classification (MUSIC) algorithms 
were investigated for spectral estimation. After the sensor nodes are deployed, they are 
grouped into clusters with a primary controlling node, the central controller, being 
established. Once the location of the sensor nodes has been determined by the central 
controller spectral estimation will be initiated and coordinated in order to evaluate the 
spectrum and determine the frequency of the SOI. An introduction into the two spectral 
estimation techniques investigated will be presented below [16, 73]. 

Spectral estimation methods estimate the frequency of the SOI impinging on our 
distributed antenna array. One of the more straightforward methods is that of the 
Periodogram. This method is widely used because the computation is efficient, as shown 
in Equation (30) [66]. In practice, a fast Fourier transform (EFT) is used to compute the 
spectrum of the given data sequence. 



(30) 


where is the Fourier transform of the received signal and is the number of 

samples used in the Fourier transform. 
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The MUSIC method was developed by Sehmidt, et a/. [47] This algorithm will be 
presented in more detail in Chapter V, Section A, because this method will also be 
investigated for determining angle-of-arrival of a signal of interest (SOI). 

The central controller will be assumed to have more processing power and a 
higher gain antenna than the individual nodes. Therefore, the central controller will be 
considered the reference antenna. The research was conducted utilizing built-in 
MATLAB® functions to compare the two methods. 

Both of these techniques have advantages and disadvantages. Research indicates 
that both perform adequately at lower frequencies, but MUSIC continues to give sharper 
results at high frequencies as seen in Figures 35 - 37. 
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Figure 35. Spectral Estimation results 157 MHz. 
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Figure 36. Spectral Estimation results 800 MHz. 
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Figure 37. Spectral Estimation results 2.4 GHz. 


As discussed earlier the periodogram takes the Fourier transform of the received 
power spectral density. As can be seen in the figures, the noise appears to have more of 
an effect at the higher frequencies. Where as in the case of MUSIC the signals are chosen 
by selecting eigenvalues that are orthogonal to the noise subspace, the effect of the noise 
has been removed. The estimates of the frequencies are taken as those values of co that 
create peaks in the spectrum. When the ensemble average of the array input covariance 
matrix is known and the noise can be considered uncorrelated and identically distributed 
between the elements, the peaks of the MUSIC spectrum are guaranteed to correspond to 
the frequencies of the signals incident on the array [16, 30, 47, 58, 62, 66, 69, 71, 73]. 

MUSIC gives better resolution for signals operating at closely spaced frequencies. 
MUSIC is highly robust, but it requires characterization of the array response. The 
periodogram has the advantage of being extremely easy to implement in FPGA hardware 
allowing for a very fast solution. However, the periodogram is unable to resolve signals 
spaced less than In! N apart [62, 66]. As described in [62], the resolution limit for the 
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MUSIC method is more eomplex, and depends on the number of samples as well as SNR 
and number of antenna elements. For more detailed deseriptions on the assumptions on 
the resolution limits of MUSIC see [62], 

However, MUSIC does not perform well in praetieal applieations where the 
number of signals present are unknown. The algorithm is unable to separate the signal 
and noise into their respeetive subspaees. Due to its simplieity and wide use it is assumed 
that the SIGINT/IW network will employ some variant of the periodogram to determine 
the frequency of the SOI [62]. 

D, ELEMENT SYNCHRONIZATION FOR THE SIGINT/IW SENSOR GRID 

An investigation was conducted into the two synchronization methods developed 
at NFS. The intent of this investigation is to explore the viability of the methods to 
synchronize the elements contained within the SIGINT/IW network. In [28], the “brute 
force” and “beam tagging” techniques were applied to a wireless radar scenario. The 
wireless Tx/Rx modules were distributed randomly over the hull of a DDx Destroyer. 

The MATLAB® code created for these simulations was altered to investigate the potential 
opportunities for use within a SIGINT/IW System. This modified code can be found in 
Appendices A and B. The random distribution of sensors in [28] was bounded by the 
dimensions of the DDx. The code was enhanced to allow for the random distribution of 
SIGINT/IW sensors over a generic random surface. The simulations in MATLAB® were 
performed to verify the “brute force” and “beamtagging” techniques. The programs were 
used to phase synchronize with 400 elements, distributed randomly. The central 
controller was located at the origin so that the transmission path length was equal to the 
norm of the x andy coordinates of element n. 

1 . “Brute Force” Synchronization Technique 

As was previously shown in Figure 9, the brute force technique functions as 
follows. All elements are initialized with zero phase and the element closest to the origin 

is selected as the reference element. Each element is selected in turn for synchronization. 

0 

When the element is selected, its phase shifter is incremented in 22.5 steps, which is 
equivalent to four-bit digitization, until the combined field is minimum. This is then 


61 



repeated for the rest of the elements. Four-bit digitization was deemed satisfactory for 
digital phase shifter quantization based on the required sidelobe level [28]. 

Figure 38 shows the phase error of each element from the reference element 
plotted against the number of iterations. A change in color denotes a new element being 
synchronized. For one realization of the 400 randomly located elements, 3247 iterations 
were required to perform synchronization. The frequency chosen was 400 MHz, which 
was the same used in [28] in order to get a good baseline comparison. 
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Figure 38. Phase Synchronization of 400 Elements using the Bruteforce Technique. 

For fixed phase steps, it is not possible to achieve complete cancellation. 

Therefore the brute force technique requires a threshold value to detect the minima when 

0 

the two signals cancel. Assuming equal amplitudes and using 22.5 steps, the final phase 
error is ±11.5 . Using phasor geometry the minimum can be determined. A threshold of 
0.2 (-14 dB) was used in [28] to determine the minima. With 16 phase steps between 0 
and 360 about eight iterations on average were required to synchronize each element. 
When the elements are synchronized with the reference element, the final phase error is 

o 

between ±11.25 as shown in Figure 39. These are the same results obtained in [28] 
showing that the algorithm is adaptable to a totally randomly distributed antenna grid. 
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Figure 39. Final Phase Error of all Elements after Phase Synehronization using the 

“brute-force technique. 


Now the frequency was changed to the three frequencies utilized throughout this 
research, 157 MHz, 800 MHz, and 2.4 GHz, to quantify the effects. As seen in Eigures 68 
-85 the frequency being synchronized did not seem to have an effect. However, these 
simulations do not account for channel characteristics, i.e., multipath, etc. that could be a 
factor, especially at high frequencies. The phase error for all cases was the same ±11.25 . 
The figures highlight the synchronization and final phase errors of individual nodes 
instead of the entire array. 
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Figure 40. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 157 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 25 square meter area. 
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Final phase error for a sample of nodes from a collection of 400 nodes in a 25 
square meter area. Frequency of the SOI is 157 MHz. 
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Figure 42. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 157 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 2500 square meter area. 


Figure 43. 
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Final phase error for a sample of nodes from a collection of 400 nodes in a 
2500 square meter area. Frequency of the SOI is 157 MHz. 
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Figure 44. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 157 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 250000 square meter area. 


Figure 45. 
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Final phase error for a sample of nodes from a collection of 400 nodes in a 
250000 square meter area. Frequency of the SOI is 157 MHz. 
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Figure 46. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 800 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 25 square meter area. 


Figure 47. 
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Final phase error for a sample of nodes from a collection of 400 nodes in a 25 
square meter area. Frequency of the SOI is 800 MHz. 
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Figure 48. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 800 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 2500 square meter area. 
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Figure 49. Final phase error for a sample of nodes from a collection of 400 nodes in a 
2500 square meter area. Frequency of the SOI is 800 MHz. 
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Figure 50. Number of iterations required for the Brute Force technique to phase 

synchronize at a frequency of 800 MHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 250000 square meter area. 
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Figure 51. Final phase error for a sample of nodes from a collection of 400 nodes in a 

250000 square meter area. Frequency of the SOI is 800 MHz. 
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Figure 52. Number of iterations required for the Brute Force technique to phase 
synchronize at a frequency of 2.4 GHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 25 square meter area. 
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Figure 53. Final phase error for a sample of nodes from a collection of 400 nodes in a 25 

square meter area. Frequency of the SOI is 2.4 GHz. 
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Figure 54. Number of iterations required for the Brute Force technique to phase 
synchronize at a frequency of 2.4 GHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 2500 square meter area. 
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Figure 55. Final phase error for a sample of nodes from a collection of 400 nodes in a 
2500 square meter area. Frequency of the SOI is 2.4 GHz. 
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Figure 56. Number of iterations required for the Brute Force technique to phase 
synchronize at a frequency of 2.4 GHz. Shown is a highlight of nodes from a 
collection of 400 nodes in a 250000 square meter area. 
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Figure 57. Phase error for a sample of nodes from a collection of 400 nodes in a 250000 

square meter area. Frequency of the SOI is 2.4 GHz. 

On average, the “brute force” technique requires half the total number of 
quantization levels to synchronize one element because the required phase shift is 
unknown. 

2, “Beam Tagging” Synchronization Technique 

The “beam tagging” method as proposed in [28] is for self-focusing or steering of 
an adaptive transmitting array. It is a technique of applying low-index phase modulation 
to one of two antennas aimed at the same target, and measuring resultant amplitude 
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modulation to correct the phase alignment between them. This technique has been used to 
phase-align lasers onto a target and for testing a large radar array [28]. The “beam 
tagging” technique can be implemented with more hardware modifications as was 
previously shown in Figure 10. As noted in [20], the key changes are the addition of a 
phase modulation circuit on the element synchronization block and an amplitude 
modulation (AM) receiver circuitry on the central controller. In each element 

synchronization block, the modem controller holds a phase shift command and is able, on 

0 

special request, to modulate the phase rapidly by ±90 from the command phase. 

The code used in [28] was again modified for the SIGINT/IW architecture. This 
code can be found in Appendix A & B. Simulation was performed to verify the “beam 
tagging” technique based on the same initial conditions used earlier. As discussed in 
reference [20], the phase corrections are performed until an opposite command is 
detected. This means that the element phase relative to reference has changed from lead 
to lag or vice versa and a balance condition is reached. The “beam tagging” technique 
detects this and terminates the synchronization cycle. 
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Figure 58. Synchronization of 400 Elements using the Beamtagging Technique. 
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The “beam tagging” technique synchronized the array after 2204 iterations, as 
shown in Figure 58, a significant reduction of 52% in the number of iterations. The 

o 

steady phase error is between ±22.5 as shown in Figure 59. The steady phase error is 
greater because when the synchronization cycle is terminated, the phase could have been 

0 o 

overcorrected by 22.5 , giving rise to the error of ±22.5 . These results are consistent with 
the results found in [28], which proves the viability of using either method to synchronize 
the individual antenna elements in a randomly distributed array. 



Element rurrbof 

Figure 59. Phase Error of all Elements after Phase Synchronization using the “beam tag” 

technique. 

In summary [28], concluded that even though the “beam tagging” technique was 
faster on the order of 2 to Sjis, the more simple “brute force” technique is preferred for 
the wireless radar. Based on the research conducted for the SIGINT/IW antenna network 
using the modified MATEAB® algorithms the same conclusions can be drawn, that the 
brute force technique will be sufficient for the SIGINT/IW sensor architecture. Therefore, 
the SIGINT/IW sensor architecture will assume that synchronization will be conducted 
utilizing some variant of the “brute force” technique. 
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This chapter provided an overview into RF signals colleetion. A basic overview 
of antenna array theory and antenna array properties was discussed. Additionally, speetral 
estimation teehniques were reviewed. Finally, the last section provided a look at the 
ongoing research by Jenn, et al. [20] into a difficult aspect of a wireless antenna. Even 
though synehronization is not central to this research it is important to give the reader an 
understanding of the complexity of this issue. Synchronization between the elements is 
assumed in this research. However, synchronization between the individual elements is a 
non-trivial requirement needed to make the SIGINT/IW wireless network a reality. 
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V. SUBSPACE METHODS 


Most direction finding antennas are configured in a uniform manner, i.e. circular, 
linear, etc. The spacing between antenna elements is critical in the design of a traditional, 
uniformly spaced phased array. This is because when the spacing exceeds certain critical 
values, grating lobes will occur [3]. Grating lobes are additional sidelobes that have the 
same amplitude as the main beam, ft is normally undesirable to have grating lobes, and in 
order to avoid them, the element spacing is required to meet the condition given by 
Equation ( 8 ). Equation ( 8 ) is rewritten as Equation (31) for convenience. 

yi , I* | . (31) 

1 + |sm 0^ I 

where d is the element spacing, is the wavelength of the highest operating frequency 
and 6*0 is the scan angle. 

Grating lobes require the coherent accumulation of the sampled signal at angles 
other than the steering angle, ft will be shown that these grating lobes and the random 
distribution of our sensors cause a correlation between the signal and the noise that make 
traditional subspace methods, i.e. MEISIC, MEM, etc. non-optimal for solving for 
direction of arrival (DOA). 

A. MULTIPLE SIGNAL CLASSIFICATION (MUSIC) 

The Multiple Signal Classification (MUSIC) algorithm introduced by Schmidt 
[47, 71] is a subspace based algorithm that exploits eigenstructure properties of the array 
correlation matrix. As discussed in references [30, 47, 69, 73], MUSIC can be used for 
estimating a range of signal parameters including the number of incident signals and their 
DOA. MUSIC is also extensively used in spectral estimation to estimate the frequency 
[73]. Earlier the spectral estimation properties of MUSIC were demonstrated. The angle 
of arrival determination properties will now be presented. Several variations of MUSIC 
exist, but in this research we will only investigate the standard form [30, 47, 69, 73]. 
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As previously stated MUSIC separates the reeeived signal plus noise into a signal 
subspaee and a noise subspace using eigenvalue decomposition [16, 30, 69, 73]. Once the 
noise subspace has been estimated, a search for angle of arrival of the desired signal is 
made by looking for steering vectors that are as orthogonal to the noise subspace as 
possible [73]. The search may be performed using either the noise subspace or the signal 
subspace over all possible DO As to determine the steering vectors. The steering vectors 
corresponding to the directional sources are orthogonal to the noise subspace and are 
found within the signal subspace [16, 30, 69, 73]. 

An example used in [71], illustrates this process. For an M element array with D 
signals incident on the array, the received input data vector may be expressed as a linear 
combination of the D incident waveforms and noise. 



b = ^>^ + n, (33) 

where 5 is the vector of incident signals, n is the noise vector, and 'F(0y) is the array 
factor or the array steering vector corresponding to the direction of arrival of the 
signal. Viewing the received vector b and the steering vector T (©y j as vectors in M 
dimensional space, it can be seen that 6 is a particular linear combination of the array 
steering vectors with through being the coefficients of the combination. Using this 
data model, the input covariance matrix can be expressed as 

=E\bb"^ = ^>E[_ss^^^^>^ +E\nn"\ (34) 

R„=E[bb‘'y'VE[ss"]'V" , (35) 

where R^^ is the signal correlation matrix £’[^55^] . 

Using eigenvalue decomposition R^g is separated into a K -dimensional signal 
subspace and a (M -7f) -dimensional noise subspace. The eigenvectors associated with 
the smaller eigenvalues of the correlation matrix correspond to the noise subspace, and 
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the eigenvectors associated with the larger eigenvalues correspond to the signal subspace. 
Estimating with a finite data sample can cause the eigenvalues corresponding to the 
noise power to not be identical. They will, however, appear as a tightly spaced cluster of 
eigenvalues with the variance decreasing as the number of data samples used to estimate 
R^g increases. Once the multiplicity of the smallest eigenvalue is determined, the number 
of signals is given by [47, 69, 71] 

D = M-K, (36) 

where M is the number of elements and K is the number of eigenvalues corresponding 
to the noise subspace. Based on the definition of eigenvalues and eigenvectors, the 
eigenvectors corresponding to the smallest eigenvalues satisfy 

for/ = Z) + l,...,M, (37) 

Observing Equation (35), we see that 

= 0, for/ = E) + l,...,M, (38) 

Since 'E is full column rank and R^g is non-singular, it follows that 

for/ = E) + l,...,M, (39) 

This implies that the column vectors of 'E are perpendicular to the eigenvectors 
corresponding to the noise. 

Summarizing the above example, the MEISIC method is based on the observation 
that the steering vectors corresponding to the DOA of the signals are orthogonal to the 
noise [16, 30, 47, 58, 62, 66, 69, 71]. Once the noise subspace as been calculated from 
the estimated array correlation matrix in Equation (35), then 'E can then be searched to 
find the eigenvectors that are orthogonal to the noise subspace. The estimates of the 
DOAs are taken as those (p that give the smallest value of 'E^ , i-e., the values 

that result in a steering vector furthest away from the noise subspace [16, 30, 47, 58, 62, 
66 , 69, 71]. The MEISIC spectrum can be expressed as 

The matrix is a projection matrix onto the noise subspace. For 

steering vectors that are orthogonal to the noise subspace, the denominator of Equation 
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(40) will become very small, and thus the peaks will occur in {(f) corresponding to 

the angle of arrival of the signal [16, 30, 47, 58, 62, 66, 69, 71]. When the ensemble 
average of the array input is known and the noise can be considered uncorrelated and 
identically distributed between the elements, the peaks of the MUSIC spectrum are 
guaranteed to correspond to the angle-of-arrival of the target signals [16, 30, 47, 58, 62, 
66, 69, 71]. However, because of the grating lobes discussed earlier the signal and noise 
are not uncorrelated between elements. Based on the earlier discussion in Chapter IV, 
Section A, Subsection I, about grating lobes it will be shown that MUSIC is non-optimal 
for the SIGINT/IW array [16, 30, 47, 58, 62, 66, 69, 71]. 

As can be seen in Figures 60-71, MUSIC is very sensitive to element spacing. 
This sensitivity will be illustrated by utilizing a variant of Equation (40) called Root 
MUSIC, as shown in Equation (41) to calculate the optimal spacing. This method is only 
applicable when using a UEA, but it provides a means to calculate the optimal spacing. 
The optimal spacing, d, for angle-of-arrival of 30° was calculated to be 0.1591549 k. In 
all of the frequencies investigated, the calculated angle-of-arrival was determined to be to 
the right of the real azimuth when the spacing was less than the optimal. Eikewise, the 
calculated angle-of-arrival was determined to be to the left of the real azimuth when the 
spacing was greater than d. The information used for the simulations is summarized in 
Tables 3-11. Tables 3-5 are for frequency 2.4 GHz. Tables 6-8 are frequency 800 MHz 
and finally Tables 9-11 are for frequency 157 MHz. 

Pmv.A'/’P^ -. <41) 

Z PiA 

i=M+\ 


Table 3. 


Number of 
Elements 

Distance 

(m) 

k = 0.125 
meters 

Erequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

4 

0.15k 

2400 

30.295° 

27.4219° 

4 

0.17k 

2400 

30.136° 

31.6406° 

4 

0.2k 

2400 

29.991° 

36.5625° 


A four element UEA with a test frequency of 2.4 GHz at three different element 


spacings. 
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Number of 
Elements 

Distance 

(m) 

L = 0.125 
meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

10 

0.15L 

2400 

30.034° 

27.4219° 

10 

0.17L 

2400 

29.993° 

31.6406° 

10 

0.2L 

2400 

29.994° 

36.5625° 


Table 4. A ten element ULA with a test frequency of 2.4 GHz at three different element 

spacings. 


Number of 
Elements 

Distance 

(m) 

L = 0.125 

meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

20 

0.15L 

2400 

29.998° 

27.4219° 

20 

0.17L 

2400 

30.018° 

31.6406° 

20 

0.2L 

2400 

29.991° 

36.5625° 


Table 5. A twenty element ULA with a test frequency of 2.4 GHz at three different 

element spacings. 


Unifoim Lintai Array wKh 4 •Itmtntt • MUSIC 



Figure 60. 2-D Array Pattern for a 4 element ULA at 2.4 GHz utilizing MUSIC. 
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Uniform Lntat Anay with 10 tfomantt ■ MUSIC 



Figure 61. 2-D Array Pattern for a 10 element ULA at 2.4 GHz utilizing MUSIC. 


Unifetm Un«st Array with 20 - MUSIC 



Figure 62. 2-D Array Pattern for a 20 element ULA at 2.4 GHz utilizing MUSIC. 
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Number of 
Elements 

Distance 

(m) 

L = 0.375 
meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

4 

0.15L 

800 

29.562° 

27.4219° 

4 

0.17L 

800 

29.936° 

30.9375° 

4 

0.2L 

800 

30.448° 

36.5625° 


Table 6. A four element ULA with a test frequency of 800 MHz at three different element 

spacings. MHz. 


Table 7. 


Number of 
Elements 

Distance 

(m) 

L = 0.375 
meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

10 

0.15L 

800 

29.975° 

27.4219° 

10 

0.17L 

800 

30.1° 

31.6406° 

10 

0.2L 

800 

30.06° 

36.5625° 

L ten element i 

ULA with a test frequency of 800 MHz at tl 

iree different 


spacings. 


Number of 
Elements 

Distance 

(m) 

L = 0.375 

meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

20 

0.15L 

800 

29.989° 

27.4219° 

20 

0.17L 

800 

30.015° 

31.6406° 

20 

0.2L 

800 

29.993° 

36.5625° 


Table 8. A twenty element ULA with a test frequency of 800 MHz at three different 

element spacings. 
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Uniform bn«ar Array with 4 elemtnt s • MUSIC 



Figure 63. 2-D Array Pattern for a 4 element ULA at 800 MHz utilizing MUSIC. 


Unilbmi Linear Ajray with 10 elements • MUSIC 



Figure 64. 2-D Array Pattern for a 10 element ULA at 800 MHz utilizing MUSIC. 
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Uniform Linear Array with 10 elements • MUSIC 



Figure 65. 2-D Array Pattern for a 20 element ULA at 800 MHz utilizing MUSIC. 


Number of 
Elements 

Distance 

(m) 

1 = 1 

meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

4 

0.15L 

157 

29.896° 

27.4219° 

4 

0.17L 

157 

29.749° 

31.6406° 

4 

0.2L 

157 

29.943° 

36.5625° 


Table 9. A four element ULA with a test frequeney of 157 MHz at three different element 

spacings. 


Table 10. 


Number of 
Elements 

Distance 

(m) 

L = 2 

meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

10 

0.15L 

157 

30.04° 

27.4219° 

10 

0.17L 

157 

30.045° 

30.9375° 

10 

0.2L 

157 

30.0024° 

36.5625° 

L ten element 

JEA with a test frequency of 157 MHz at t 

u’ee different 


spaeings. 
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Number of 
Elements 

Distance 

(m) 

X = 2 

meters 

Frequency 

(MHz) 

DOA 

RootMUSIC 

Azimuth 

Calculated 

20 

0.15L 

157 

30.004° 

27.4219° 

20 

0.17L 

157 

29.968° 

30.9375° 

20 

0.2>. 

157 

29.992° 

36.5625° 


Table 11. A twenty element ULA with a test frequency of 157 MHz at three different 

element spacings. 


Unfoim ljn*ai A/ray wtlt 4 alanxnlt . MIJSIC 



Figure 66. 2-D Array Pattern for a 4 element ULA at 157 MHz utilizing MUSIC. 
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Unifomn Linear Array with 10 elements • MUSIC 



Figure 67. 2-D Array Pattern for a 10 element ULA at 157 MFlz utilizing MUSIC. 


Uniform Linear An?y vnlh 20 elements - MUSIC 



Figure 68. 2-D Array Pattern for a 20 element ULA at 157 MHz utilizing MUSIC. 


The results from the above simulations were taken at frequencies of 157 MHz, 
800 MHz and 2.4 GHz. The antenna element spacings were normalized based on 
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wavelength. The largest element spaeing distance was 0.4 meters apart for the 157 MHz 
frequency. The closest spacing was 0.01875 meters for the 2.4 GHz frequency. These 
spacings are very unrealistic for an antenna field that is expected to be deployed from 
artillery shells, dropped from UAVs or scattered via a UUV. Now a look at the array 
pattern characteristics when the elements are spaced at wavelengths equal to or greater 
than half-wavelength were investigated. The data for these simulations can be seen in 
Tables 12 - 14. Figures 69-71 illustrate that the inconsistency in the array pattern is 
much worse as the element spacings increase. 


Number of 
Elements 

Distance 

Frequency 

40 

;l 

2 

2.4 GHz 

40 

3A 

2 

2.4 GHz 

40 

2A 

2.4 GHz 


Table 12. A forty element ULA with a test frequency of 2.4 GHz at three different element 

spacings equal to or greater than half-wavelength. 


40 element ULA 


UndormltmNK Aiinywnh lOetnnwnli ML18IC 
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Figure 69. 2-D Array Pattern for a 40 element ULA at 2.4 GHz utilizing MUSIC where 

the elements are spaced at half-wavelength or greater. 
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Number of 
Elements 

Distance 

Frequency 

40 

;l 

2 

800 MHz 

40 

31 

2 

800 MHz 

40 

21 

800 MHz 


Table 13. A forty element ULA with a test frequeney of 800 MHz at three different element 

spacings equal to or greater than half-wavelength. 


40 element ULA 
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Figure 70. 2-D Array Pattern for a 40 element ULA at 800 MHz utilizing MUSIC where 

the elements are spaeed at half-wavelength or greater. 


Number of 
Elements 

Distance 

Erequency 

40 

1 

2 

157 MHz 

40 

31 

2 

157 MHz 

40 

21 

157 MHz 


Table 14. A forty element ULA with a test frequency of 157 MHz at three different element 

spacings equal to or greater than half-wavelength. 
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Figure 71. 2-D Array Pattern for a 40 element ULA at 157 MFlz utilizing MUSIC where 

the elements are spaced at half-wavelength or greater. 

When more realistic distances are used it can be seen in the above figures how the 
results become very erratic. Other variants and/or adaptations of MUSIC may be able to 
better cope with this phenomena but these will not be explored in this research. Now 
we’ll take a look into the Maximum Likelihood method. 

B. MAXIMUM LIKELIHOOD 

The method of Maximum Likelihood (ML) estimation is one of the most popular 
estimation methods used in signal processing. For a DOA estimation problem, the term 
“maximum likelihood” is used for the method that finds the ML estimate of the direction 
rather than of the power [16, 33]. This method uses the array weights, which are obtained 
by minimizing the mean output power in the look direction [16, 30, 47, 58, 62, 66, 69, 
71]. 

An expression for the power spectrum is given by 

Just like with MUSIC due to the grating lobes discussed earlier in Chapter IV, 
Section A, Subsection 1, MLM is non-optimal for the SIGINT/IW array. This can be 
easily seen in Figures 72-83. 
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Unifoim bneat Array with 4 tlemenls - MAXIMUM LA<EUHOOD METHOD 



Figure 72. 2-D Array Pattern for a 4 element ULA at 2.4 GHz utilizing MLM. 

Uniform Linear Array with 10 akmenu - MAXIMUM LIKEUHOOD METHOD 



Figure 73. 2-D Array Pattern for a 10 element ULA at 2.4 GHz utilizing MLM. 
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IJoiform Uneat Airay with 20 «l«m«nt« - MA30MUM UKEUHOOO METHOD 



Figure 74. 2-D Array Pattern for a 20 element ULA at 2.4 GHz utilizing MLM. 


Uniform Linear Array vnth 4 elements - MAXIMUM UKEUHOOO METHOD 



Figure 75. 2-D Array Pattern for a 4 element ULA at 800 MHz utilizing MLM. 
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Uniform Lii>**r Air«y vnth 10 titmontt • MAXIMUM LIKELIHOOD METHOD 



azimutli anjl^-of arrival ^ in dagraas 


Figure 76. 2-D Array Pattern for a 10 element ULA at 800 MHz utilizing MLM. 


Unrfam Unear Array vuith 30 elements - MAXIMUM UKEUHOOO METHOD 



Figure 77. 2-D Array Pattern for a 20 element ULA at 800 MHz utilizing MLM. 


As in the case for MUSIC, the results from these simulations were taken at 
frequencies of 157 MHz, 800 MHz and 2.4 GHz. The antenna elements were spaced at 
fractions of wavelengths apart for convenience of computation. The largest spacing was 
0.4 meters apart for the 157 MHz frequency. The closest spacing was 0.01875 meters for 
the 2.4 GHz frequency. 
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Unrterm Linear A/ray vnth 4 elements - MAXIMUM UKEUHOOC METHOD 



Figure 78. 2-D Array Pattern for a 4 element ULA at 157 MHz utilizing MLM. 


Uniform Linear Anay verth 10 elements • MAXIMUM LIKELIHOOD METHOD 



Figure 79. 2-D Array Pattern for a 10 element ULA at 157 MHz utilizing MLM. 
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Uniform Linear Array wrth 20 elements - MAXIMUM UKEUHOOD METHOD 



azimuth angle-of-arrival <{> in degrees 


Figure 80. 2-D Array Pattern for a 20 element ULA at 157 MHz utilizing MLM. 



Figure 81. 2-D Array Pattern for a 40 element ULA at 2.4 GHz utilizing MLM. 
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Figure 82. 2-D Array Pattern for a 40 element ULA at 800 MHz utilizing MLM. 


Figure 83. 


40 element ULA 
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2-D Array Pattern for a 40 element ULA at 157 MHz utilizing MLM. 


As in the ease of MUSIC, when the elements become further apart the results 
become more erratic. Other variants and/or adaptations of MLM may be able to better 
cope with this phenomena but these will not be explored in this research. Therefore, it 
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will be impractical to impossible to use traditional subspace methods to determine the 
angle of arrival. Hence, this research will investigate methods that can mitigate the 
grating lobe problem. 

This chapter has provided insight into the problems with using traditional 
subspace methods with the SIGINT/IW wireless sensor array. Since the sensor nodes will 
in some cases be spaced greater than a half-wavelength apart, the grating lobes cause the 
signal and noise to become correlated rendering the traditional subspace methods 
inadequate for determining the direction of arrival. Therefore this leads into the next 
phase of the research, which is to find a method suitable for determining this direction of 
arrival in order to utilize beamforming to enhance the collection against the target. 
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VI. DIRECTION-OF-ARRIVAL 


This chapter provides a baekground into the two methods analysed in this 
dissertation for line of bearing determination, time differenee of arrival and 
beamforming. Both methods were analyzed in order to determine the most expeditious 
and energy effieient method to determine a line of bearing to the target. This ehapter will 
give the reader an appreeiation into some of the eonstraints and tradeoffs assoeiated with 
both methods. 

A, TIME DISTANCE-OF-ARRIVAL (TDOA) 

Time difference of arrival (TDOA) takes advantage of the fact that a transmitted 
signal will arrive at the different elements at different time instants. It will be shown that 
TDOA is not affeeted by the grating lobes, whieh bemoans the subspace methods of 
MUSIC and MLM. In faet it will be shown that as the baseline d , between the sensor 
node and the eentral eontroller inereases the resulting TDOA measurement varianee gets 
smaller and subsequently the resulting TDOA gets better. In the SIGINT/IW wireless 
antenna array eaeh sensor forms a two-node eolleetor pair with the central eontroller, as 
shown in Figure 84. The basic principle of the two-element linear array will then be 
extended to eneompass the entire antenna grid. The sensors within our grid will take 
advantage of the faet that our grid is densely populated and that the sensors are spread 
over a wide area [29]. 

Defining TDOA in terms of the emitter and antenna element positions is 
aeeomplished using the vectors n and from Figure 84 where £ is the distanee from 
the eentral eontroller to the emitter and is the distanee from an individual antenna 
element to the emitter [29]. 
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Figure 84. 


y 



TDOA Geometry for a eolleetor pair. [After Ref 47] 


Equation (44), illustrates the differenee in the length of the two veetors 
2 J. = [xg - / 2; and ^ J / 2; y , divided by the speed of light c, yielding the 

TDOA; 


TDOA = 



1 

c 



+ ye 



(43) 


where x^ and are the respeetive x-y eoordinates of the emitter and d is the distanee 
between the two-node eolleetor pair eommonly referred to as the baseline. 

When the two-node eolleetor pair are in motion, multiple TDOA measurements 
ean be used to estimate emitter position in a manner analogous to a navigator performing 
triangulation. For the wireless antenna array network the two-node eolleetor pairs are not 
in motion therefore only one measurement was taken from eaeh two-node eolleetor pair. 

The narrowband signal shown in Equation (Q-10) will be assumed to be the 
signal-of-interest arriving at our two different sensors. The likelihood funetion for the 
reeeived signal is shown in Equation (44) and is derived in detail in Appendix Q [72]. 


100 



E[?{t)\ = AI, 


^2Aq'^ 
V^o y 


exp 


-1 


2N, 


J|r(^)| +A^s^^^{t)dt 


0 T 
2 


q = 




( 44 ) 


where is the amplitude of the reeeived signal, f (t) is the reeeived signal’s eomplex 


envelope plus noise as defined in Equation (Q-16). (t) is the eomplex amplitude of the 

reeeived signal. is the noise power speetral density. 


The signal is assumed to arrive in an interval between -T12 and T12. /(,(□) is 

the modified Bessel funetion order zero. This funetion is a result of faetoring out the 
earner phase as a nuisanee parameter. The earrier phase is assumed to be a uniform 
random variable between -n and n. An assumption is made that the SNR is high. The 
argument of the exponential is assumed to be eonstant and the only thing assumed to be 
unknown is the argument of the modified Bessel funetion. Therefore the likelihood 
funetion is maximized when the argument of the modified Bessel funetion, speoifieally, 
q is maximized [72]. 

It is assumed that the envelope of the reeeived signal is a elose approximation to 
the transmitted signal. Therefore, q ean be expressed as 


q = 


J {t-T^dt 


(45) 


where the reeeived signal has been rewritten in terms of the transmitted signal from 
Equation (6) only delayed in time, r . As shown below 




( 46 ) 
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The high frequency components have been removed during the filtering process. 

It is assumed that a quadrature detector with a matched filter is used. The time at which 
the envelope peaks is the approximate value of the delay [72]. 

The Cramer-Rao bound for the variance of this estimation taken from [72] will 
now be derived for the narrowband signal case. The received signal is assumed to be 

r{t^ = As^{t-T^ + h{t^, (47) 

The derivation uses the large signal approximation. Now taking the natural log of 
likelihood function results in the following; 


ln[j!7(r|r) =ln(y4) + ln 
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+ ln 

exp 
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\ ^0 J 

/ 





( L 
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2N, 




0 r 
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, (48) 




when the argument of the modified Bessel function is maximum then 


In/g (x) = X . 

Based on this assumption Equation (48) can be rewritten in the form 


(49) 


ln[;7(r|r)]^ln(^) + ^ + 


T 

^ \ + AWt)dt 

^^^0 T_ 

This term is approximately independent of t. 


(50) 


Due to its independence on r , the last term in Equation (51) can be removed. 
Substituting back in for q from Equation (47), Equation (51) now becomes 


ln[/»(r|r) =ln(T) 


+ - 


A 


A 

I r{t)s* {t-T)dt 




I 

I {t- + h{tY^sJ [t- dt 


( 51 ) 
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As discussed in [72], because of the large SNR assumption the noise term is 
dropped and r » 0 in the last term. Therefore Equation (51) is rewritten as 

T 

ln[j!7(f |r)]^ln(^) + ^ J As,[t -t)s* [t)dt , (52) 

2 

An important aspect of calculating the time difference of arrival is the variance 
associated between multiple measurements. This variance determines the plus or minus 
error associated with the subsequent line of bearing determination. Equation (52) 
provides convenient representation for solving for the Cramer-Rao bound of the 
associated time difference measurement variance. Now a more detailed discussion of the 
cross-ambiguity function will be offered. The cross-ambiguity function will be utilized to 
calculate the time difference of arrival. 

1, Cross-Ambiguity Function (CAF) 

The complex envelope of the transmitted signal is collected by the sensors then 
the sampled signals are sent back to the central controller where the cross-ambiguity 
function (CAF) is used to determine the TDOAs. An FDOA, can also be determined from 
calculating the CAF, which could be very useful if prosecuting a mobile target [24]. 
However, our target is assumed stationary. The complex envelopes of the two signals will 
be denoted by s, {t - r) and s* (t), where t is the time of the first sample, and r is the 
time lag between the two. The CAF is defined as 

T 

CAF{t,w) = ^ \ 

T 

CJF(0) = l = ^|i,(()y(r)d( , (53) 

Iht ^ 

2 

T 

1 i 2 

= 2 J" 
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where is the energy in the transmitted signal, / is frequency, T is the integration 

time, and * denotes the complex conjugate. The TDOA and FDOA can be calculated 
simultaneously, since the function peaks when r = TDOA and / = FDOA [56]. The 

solution is determined by locating the peak of the CAT, |C4F (/, r)| [24, 56]. It is worth 

noting that when the sensors or the emitter are not in motion that the cross-ambiguity 
function is simply the autocorrelation function. In Figures 85 - 88, the TDOAs, of the 
sensor located the farthest from the central controller and the sensor located closest to the 
central controller are depicted. The frequency difference of arrival is zero which is 
expected. 

Cioes Anibiguity Function 

4 

X 10 



Figure 85. 3-D plot of the CAF at the maximum distance from the reference node. 
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Figure 86. 



2-D plot of the CAF at the maximum distance from the reference node. 
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Figure 87. 3-D plot of the CAF at the minimum distance from the reference node. 
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Figure 88. 2-D plot of the CAF at the minimum distance from the reference node. 

The CAF is computationally expensive. Three methods were explored in [23] for 
computing the CAF. Due to the concern for processing capability as it pertains to power 
management the more efficient algorithm discussed in [23] was utilized. This algorithm 
was designed for use with mobile collectors and was adapted for use with the wireless 
SIGINT/IW antenna network. The modified MATLAB® code is in Appendices C - D 
[29, 32, 42]. 

Equation (52) can be rewritten in terms of the cross-ambiguity function as follows 

'lA^F 

ln[;7(F|r)]-ln(^) + —^|C^F(r,0)| , (54) 

where co is zero since it is assumed that both the sensors and the target are stationary or 
slow moving. As defined in [72] the Cramer-Rao bound has been defined for time 
estimation as 
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dHn[p[y\x)\ 

dx^ 


(55) 


For a more complete derivation of Equation (55) see Appendix T . Rewriting 
Equation (55) in terms of the CAP , the Cramer-Rao Bound becomes 


Var [f - r > 


-1 



CAF(t) " 

I2j 

df 


(56) 


In [72] the second partial derivative of the CAP is shown in Equation (57) 


a' CAF[t) 


dP 


-1 




n2 


J(2^/)'|^,(/)f #+ \{2nf)\sXf)^ df 




(57) 


It can be clearly seen how the baseline distance affects the variance of the measurements. 
As the baseline distance increases so does the result of the CAP as depicted in figures 
113 - 116. Since the CAP is in the denominator of the variance as seen in Equation (57), 
as the CAP gets larger the variance decreases and vice-versa. 

A radian bandwidth term was defined in [72] to simplify Equation (57) 


Far[f-r > 


-1 


2A^P. 


(58) 


‘ /3^ 




where pd is the radian bandwidth and is defined as follows: 




— J I {2nf)^ \s, (/)f df - j^— J J [2nf)\s, [ff df , 


(59) 


If 1^, (/)| is assumed to be a rectangular spectrum then equation (59) can be 
rewritten as shown in Equation (60) 
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where the speetrum is defined -H B^< f <\l B^. After the integration reduces to 




3£. 


After substitution the Cramer-Rao bound becomes 


( 61 ) 


cr 
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2^; 

V y 


(62) 


Taking the square root and then canceling like terms 


^TDOA 




where = —, 

° w 


(63) 


\ \ N ) 

N is defined as the noise power and W is the bandwidth of the receiver. The denominator 
will be further refined as follows: 


v3 ^ 

^TDOA ~ ^ where; 7 =——, (64) 

where ?; is a processing gain signal-to-noise ratio taken at the output of the receiver that 
is based on the signal-to-noise ratio of the received signal, integration period and the 
bandwidth of the received signal. B^ is the bandwidth of the signal-of-interest. Therefore 
the Cramer-Rao bound can be approximated 

^ TDOA ~ „ • (65) 

B,ri 

For multiple estimates the variance becomes 
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( 66 ) 


Bsl j 

where Y is the number of estimates taken. For our purposes this will be the number of 
sensors that participate in the solution. The importance of this variance will be discussed 
in more detail in the next section. 

2. Newton Raphson Technique 

Once the TDOAs were determined from the two-node collector pairs formed by 
the central controller and the sensor nodes the Newton-Raphson technique was used to 
calculate a geolocation position. The Newton-Raphson algorithm is based on estimation 
theory and uses an over-determined set of linear equations of the form 


m =f(^ t), (67) 

where m is a k length vector representing the TDOA measurements and z is the 

location of the emitter. In the case of the SIGINT/IW sensor network m will contain the 
TDOAs calculated from each of the two-node collector pairs and the central controller. 
The size of the vector will be equal to the number of sensors, Y, used to form the two- 
node collector pairs. The number of measurements is often much greater than the number 
of position terms. 


It is envisioned that the time difference of arrival will be calculated using the 
cross-ambiguity function discussed earlier. The TDOA measurements in the m vector are 
compared with the TDOA estimates in the m. vector, where the m. vector is the TDOA 
measurements for the current estimated position as shown in Equation (67). 

Sm = m-m.= = AS^^^ (68) 

By iterating this process the current solution is used as the next estimated emitter location 
until the solution falls below some predetermined small threshold. 

z-=^A‘WA~^ ^A‘W^Sm , (69) 

A least squares estimation is applied to yield an emitter position estimate. 


no 



3, Timing Variance 

Now a discussion on the uncertainty in the TDOA estimates will be offered. 
Uncertainty in each TDOA measurement is caused by measurement noise and translates 
into uncertainty in isochrone position creating an isochrone width. The isochrone width 
grows with distance from the central controller. The appropriate branch depends upon the 
sign of the TDOA [29, 32, 42, 67]. 

The uncertainty is expressed as a measurement variance in Equation (65). The 
variance, may increase due to ionospheric, atmospheric, or other effects. As will 

be shown later in this chapter the TDOAs are very sensitive to increases in this variation. 

As can be seen in Figure 89, a sufficient receiver bandwidth is required in order to 
detect signals at a realistic SNR. It will be shown that the red curve (top), which 
corresponds to a measurement standard deviation of 1 ns is as low as the measurement 
standard deviation should be in order to maintain the accuracy depicted above. 
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When referring to Equation (65) and Figure 89, it can be seen that the only 
parameters that a designer can control is the receiver bandwidth and the noise figure of 
the receiver, which minimizes the amount of noise the receiver adds to the received 
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signal. The signal power of the reeeived signal is a funetion of many variables, i.e., 
distance from the source, multipath channels, atmospherics, etc. In all cases these are 
outside the control of the designer. Therefore the receiver bandwidth and noise figure of 
the receiver have to be chosen in order to provide a measurement standard deviation at or 
below 1 ns. 

As a point of reference, signal-to-noise ratios have been calculated utilizing the 
Cramer-Rao bound Equation (65) for some COTS receivers, as seen in Table 15. 


Model 

Frequency 
Range (MHz) 

<JtDOAv(s) 

Receiver 3 dB 
Bandwidth 
(MHz) 

Signal-to- 
Noise Ratio 
(dB) 

AD 8347 Direct 
Conversion 
Quadrature 
Demodulator 

800 - 2700 

le-9 

90 

15.7 

AD 8348 
Quadrature 
Demodulator 

50- 1000 

le-9 

500 

0.8278 

DRT 4011 
Wideband Tuner 

10-3000 

le-9 

30 

25.3 

Teams entinel 
Signal 
Acquisition 
Erontend 
(V/UHF) 

30-300 

le-9 

40 

22.8 


Table 15. Required received signal-to-noise rations (SNR) for a few commercially available 

receivers. 


The AD 8347 Direct Conversion Quadrature Demodulator provides coverage 
within the frequency range utilized by more modern communications. The AD 8348 
Quadrature Demodulator would give a better coverage range for maritime and classic 
military communications, i.e., push-to-talk, etc. 

The measurement noise covariance is also affected by the number of TDOA 

measurements utilized, which can be seen in Equation (66). The accuracy can be 

increased by the number of TDOAs, estimates per two-node collector pair, Y that are 

incorporated into the line of bearing estimate. As stated in [28] the accuracy of the line of 

bearing for a two-dimensional case estimate is as follows: 
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( 70 ) 



covariance becomes smaller the less effeet it has on the aeeuraey of the line of bearing 
estimate. Likewise as the number of estimates per two-node eolleetor pair, are inereased 
the less impaet they have on redueing the eovarianee. The question then beeomes how 
many TDOA, two-node eolleetor pairs, estimates are needed and whieh sensors should be 
seleeted, because if we only need so many TDOA estimates the eentral eontroller needs 
to select the best ones. There are several methods of seleeting the TDOA estimates or 
whieh sensors should be utilized. Seleetion methodologies is a signifieant researeh area in 
determining geoloeation and more information ean be found in [27] and [32]. For the 
purposes of this researeh the 15 estimates per two-node eolleetor pair will be used. 


The previous discussion is important to understand beeause the least number of 
estimates required in formulating a good solution means the less transmissions required 
aeross the network and less proeessing for the sensors and eentral eontroller, which 
directly translates into power and bandwidth savings. It ean be surmised that the more 
data available to the algorithm, the more aceurate the estimate as shown in the above 
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figures. The fewer the number of TDOA estimates, the less aeeurate the line of bearing 
estimate. The higher the number of TDOA estimates, the more aeeurate the line of 
bearing estimate. The tradeoff is that ealeulating the TDOAs is a computationally 
expensive process, so the greater the accuracy required the more TDOAs that have to be 
calculated. This relationship can be shown that for a given timing variance, (J^tdoa > 


Y estimates are required to generate the TDOA estimates. Beyond the Y TDOA 
estimates the decrease in diminishes to the point where no real advantage is 

gained. The relationship is as follows: 

^ TDOA ’ (' 74 ) 

where Y is the number of TDOA estimates, (J^j-doa the timing variance. Three or more 


TDOA estimates are all that is required for the Newton-Raphson Algorithm to compute a 
solution. Since the accuracy is inversely proportional to the covariance and the 
covariance goes to zero as 1/Y, then the more estimates that are used, the more accurate 
the solution. However, the amount of increase for each additional TDOA estimate would 
decrease and the noise measurement variance would become smaller and smaller as the 
number of estimates went up. This can be seen clearly in Figures 90-91 below. 
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Noise Variance vs. Number of Estimates 


o 



(b) 


Figure 90. Measurement Noise Variance, (a) Receiver Bandwidth 30 MHz. (b) Receiver 

Bandwidth 40 MHz. 


116 







































































(a) 


117 








































































Measurement Noise Variance vs. Number of Sensors 


o 


Estimates, Y 
(b) 

Figure 91. Measurement Noise Variance, (a) Receiver Bandwidth 90 MHz. (b) Receiver 

Bandwidth 500 MHz. 



Power Ratio 

dB 

20 

13.01 

40 

16.02 

50 

16.99 

70 

18.5 

80 

19.03 

100 

20 


Table 16. Received signal-to-noise power ratio conversion to decibels. 


It can be seen from figures 90-91 how the received SNR lowers the number of 
estimates required for a given noise measurement variance. Table 16 provides a 
convenient conversion between power ratio and decibels. Likewise it can be seen from 
the figures how the received bandwidth also affects the number of estimates required. It 
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is interesting to note that the higher the reeeiver bandwidth then the number of nodes, 
TDOA estimates, has less of an impact. However, in all cases it can be seen that above 10 
-15 estimates the decrease in noise measurement variance begins to wane. The accuracy 
benefit achieved for an increase above 10-15 estimates is too computationally 
burdensome for the amount of decrease in the measurement variance. As previously 
mentioned the TDOA estimates are computationally expensive to calculate and 
minimizing this cost will assist in speeding up the solution of a line of bearing as well as 
reducing the processing load on the central controller. For example, using the t ic-toc 
command in MATLAB® it was determined that it required 2 secs/two-node collector pair 
to calculate TDOA using the cross-ambiguity function. Using only 15 estimates requires 
30 seconds whereas using 80 estimates would take up to 160 seconds, therefore using 
only 30 estimates, a 81% faster solution is provided at a noise variance loss of 

3(10) (at 30 MHz bandwidth). Even fewer sensors could be used but as previously 
stated the accuracy of the line of bearing suffers. 

The fact that the grid is composed of several spatially spread two-node collectors 
will be leveraged to resolve the ambiguous bearing and calculate a line-of-bearing to 
begin beamforming. As stated earlier the CAF is computationally expensive. Even 
though the computations would take place at the central controller the individual nodes 
would need to transmit the requisite data to the central controller. 

4, Line-of-bearing Estimation 

The selection using this method is described as follows. The central controller 
will randomly select Q nodes that are located the furthest from the central controller. 

The central controller will then form a two-node collector pair with each of the selected 
nodes. 

The nodes will begin transmitting the received data back to the central controller. 
The central controller using the CAE, will calculate the time difference of arrival. Eor this 
method the nodes will send multiple received data samples back to the central controller. 
The central controller will then calculate multiple TDOA estimates for each two-node 
collector pair. 
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The total transmissions are defined as 

TrooA=^^, ( 75 ) 

where fl is the total number of two-node eolleetors and Y is the total number of 
reeeived data samples sent baek to the eentral eontroller per two-node eolleetor pair. 
These reeeived data samples are then used by the Central Controller to determine the line 
of bearing. 

If the cost per transmission is S joules / transmission , X\ joules is the energy 
expended when the receiver is on, but waiting to transmit and assumes the worst case. 

X 2 joules is the energy required to fully energize a transmitter that has been in a sleep 

mode and also assumes the worst case. The total energy cost to calculate a line of bearing 
is then defined as 

^TDOA ~ ^'^TDOA “ l) ^X\ + ^Xl 

= ^Y + (Q-1)Y(0.1^) + Q(0.5^)’ 

where is in units of joules. The first term in Equation (76) is the energy expend the 

data. The second term in Equation (76), is the total energy expended waiting to transmit. 
This term assumes that a node awaits for all other nodes to transmit before it can 
transmit. X\ is defined in terms of a fraction of the energy expended to transmit the data. 
The third term is the energy expended to bring all of the participating nodes out of a sleep 
mode, which assumes that all participating nodes were in a sleep mode. Xi is also 
defined as a fraction of the energy to transmit the data. The energy cost per sensor is 

S,„,=<Sr + Y(0.1<J)-E^ + (0.5<J), (77) 

where is measured in units of joules/sensor. 

A MATEAB® algorithm developed by Eoomis, et al. at the Naval Postgraduate 
School was utilized. The MATEAB® routine was designed for use with high-speed 
mobile collectors and had to be modified for use with the wireless SIGINT/IW antenna 
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network. The algorithm is designed to geolocate emitters based on the TDOA 
measurements. The modified code can be found in Appendices F - G. 

Now a specific example will be shown. The data for this example is shown in 

o 

Table 17. The target is at a 45 bearing. As can be seen in Figures 92 - 94 the accuracy of 
the line-of-bearing estimate increases as the number of estimates from each sensor 
increases, which can be seen by the convergence of the line-of-bearings towards the 
target. This is consistent with the results from subsection 3 above that as the number of 
estimates increases the variance will decrease by a factor of 1/Y. 


Freq 

(MHZ) 

# of Trans,, 

Ttdoa 

Area, A 2 

(m') 

Total 

Energy, 

Etdoa 

(joules) 

# of TDOA 

estimates, 

Ti 

Energy per 

sensor, H 

(joules/sensor) 

157 

40 

2500 

455 

10 

11.255 

157 

60 

2500 

66.55 

15 

16.6255 

Tab] 

e 17. Multiple TDOA Estimates. 


Lne of Bearing 



«position in meiois 


Figure 92. Isochrones from 4 two-node collectors with 5 TDOAs/two-node collector pair 

against an emitter at 45“ 
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Ui>e of Bearing 



Figure 93. Isochrones from 4 two-node collectors with 10 TDOAs/two-node collector 

pair against an emitter at 45° 


bne of Bearing 



Figure 94. Isochrones from a 4 two-node collectors with 15 TDOAs/two-node collector 

pair against an emitter at 45° 


It can be observed in the above figures that the isochrones from the two-node 
collector pairs begin to converge towards the target. As expected as the number of 
estimates increases the lines-of-bearing get more accurate. 
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This method would obviously be able to provide a line of bearing to the target to 
facilitate follow-on beamforming for enhanced collection. Consequently, a proposed 
approach for a wireless antenna sensor grid using TDOA to determine the line of bearing 
to the target is as follows: 

For this method the individual steps are listed below: 

Step 1. The central controller acts as a reference antenna by gaining initial 
intercept and frequency-of-interest (FOI) determination. 

Step 2. The central controller will then coordinate the direction-of-arrival (DOA) 
determination using n number of sensors. The central controller will calculate Y TDOA 
estimates per fl two-node collector pairs. The central controller will utilize the cross¬ 
ambiguity function presented earlier to calculate the TDOAs between each of these Q 
sensors and the central controller. 

Step 3. The central controller will then use the modified Newton-Raphson 
technique to determine a line-of-bearing to the target. 

Step 4. The central controller will then begin the beamforming process. Section B 
will cover in detail two various methods for adaptive beamforming. 

The TDOA determination process would take a total of seconds. 

tiDOA ~ C2s + OYtj + 0^2 5 

where Q is the number of nodes, ^ is the cost per TDOA estimate to perform the CAF, 
Y is the number of TDOA estimates per two-node collector pair and s is the time 
required to calculate the line of bearing using the NRT. The time a node takes to transmit 
is tj. This term assumes a TDMA type medium access, in which a node has a time slot 
and awaits for all other nodes to transmit before it can transmit. The time required to 
awake the nodes is and assumes the worse case, that all nodes selected are in a sleep 
mode. 

A method for computing TDOAs and then subsequent lines-of-bearing has been 

presented. As was shown the CAF and Newton-Raphson method as implemented was 

able to determine the line of bearing to the target. Now the central controller would be 
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able to adaptively form beams in the direetion of the target to enhance collection and/or 
offensive information operations. In the next section, beamforming will be considered as 
an alternative for discussion. 

B, ADAPTIVE BEAMFORMING 

The primary intent of using beamforming is to enhance collection and/or 
offensive information operations by taking advantage of the number of elements and 
having them work coherently to increase their respective detection and transmission 
ranges. However, as mentioned earlier the azimuth of the target needs to be known a 
priori in order for the central controller to know where to point the beam. In Section A, a 
proposal was offered for using TDOA to determine the target azimuth. Now a discussion 
will be presented on using beamforming as a method of determining the azimuth of the 
target by scanning the beam analogous to a radar searching for the azimuth with the 
highest received power. 

The beamformer method is the most basic method employed in direction of 
arrival estimation with an array antenna. A central controller can scan the beam 
electrically in a 360° arc, wherein the main beam scans and determines the direction in 
which the output power of the array becomes a maximum [18, 69]. The directions that 
give large power outputs can then be taken as the estimated direction of the target signals. 
When the power is plotted against these directions, it exhibits a peak in the direction of a 
target signal. Depending on the type of beamformer used, many different methods can be 
used [69, 71]. 

Almost all DOA estimation algorithms, whether narrowband or wideband, use 
narrowband beamforming techniques to estimate the DOA for the many different 
frequency bands. These separate estimates are then combined to get one estimate based 
on feasible statistical observations. [71] An example to illustrate this point will build on 
the discussions from Chapter IV, Section B, subsection 3. The narrowband beamformer 
from Equation (22) is rewritten for convenience and becomes Equation (79). 

= ( 79 ) 
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where the weighting veetor w emphasizes one particular direction. Given samples 
h(l),h(2),...,h(A'^), the output covariance is 

1 V T 1 V 

tv t=i tv 

Where is obtained as 

1 ^ 

( 81 ) 

tv t=i 

Many different choices of the weighting vector w can be made leading to different 
properties of the beamforming schemes [62, 70]. 

As discussed in [5] and [62] the weighting vector for the conventional 
beamformer is chosen in order to maximize the received power in a certain direction ^ , 
as 


'?(«») 

Wrf = , 

the classical spatial spectrum is obtained 


P.fW 


'v''WRss'pW 

Y“ (MW 


This is commonly referred to as the Bartlett beamformer [5,58, 62]. 


(82) 


(83) 


As discussed earlier, our wireless phased array will consist of SIGINT sensor 
nodes that are deployed randomly over an area of interest, resulting in random 
positioning of the sensor nodes. Subsequently, these sensor nodes form clusters and their 
locations are determined using location discovery techniques. These locations are then 
reported back to the central controller. In addition, the SIGINT/IW antenna elements will 
in some cases not be located within a half-wavelength from each other. Therefore an 
algorithm capable of mitigating the grating lobes discussed in Section A, subsection 1 
will have to be employed. The array will need to be able to scan 360 searching for target 
signals. Because of these factors the pattern of the array will need to be controlled 
dynamically. Therefore, an adaptive algorithm will be implemented. The azimuth that 
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gives the largest power output ean then be taken as the estimated DOA of the ineident 
signals [66, 69]. 

Once the target has been localized the central controller can enhance the gain 
toward the desired signals and can null toward undesired signals. Correct control of 
beams and nulls is indispensable for the adaptive antenna. The pattern of the array will be 
controlled by dynamically varying the phase and amplitude of the received signal of each 
element. 

Various adaptive algorithms exist, such as the least mean squares (LMS), 
recursive least squares (RLS) and constant modulus algorithm (CMA). The LMS 
algorithm recursively computes and updates the complex weights. [27, 73] The LMS 
algorithm requires knowledge of the desired signal, which will be assumed. This research 
will build on the research conducted in [68] using the LMS algorithm. Figure 95 shows 
the block diagram of a narrow-band adaptive beamformer using the LMS algorithm. In 
addition to being able to steer the main beam towards the signal of interest, the LMS 
algorithm can also simultaneously suppress interfering signals through the calculation of 
the complex weights. The grating lobes discussed earlier are also a consideration in 
determining which method to employ in determining the DOA. It will be shown that the 
LMS algorithm overcomes the effect of these lobes by eliminating all periodicities in the 
element locations by randomly selecting nodes to form the beam. In [57] this process is 
referred to as thinning the array. 
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Figure 95. MxN narrow-band adaptive beamformer. [From Ref. 27] 

The weights are ealeulated by analysis of the error between the observed and 
desired signal. These weights are used to adaptively form the beam in the desired 
direetion [69, 73]. The error is defined to be: 

e(t) = a(t)-6(t), (84) 

where a(t) is the referenee signal and 6(t) is the output of the beamformer. In praetice, 
a(t) needs to be elosely eorrelated to the desired signal 5 (t). 

As the signal is reeeived by the individual array elements the eomplex weights are 
iteratively adjusted and applied. The eomplex weights are ehosen to minimize the mean- 
square error between the beamformer output 6(t) and the referenee signal for eaeh 
iteration. This is shown in Equation (85) [27, 69, 73]. 

E[e\t)\ = E [a{t)-b{t)f , (85) 

whieh can be expanded to the form given by 

E^e^ = e\jx^ -2w^E'[«(t)5(t)] - l - E\^s{t^s^ (0]lf ’ 

where E [□] is the expectation operator. The minimum mean-square error (MMSE) is 
obtained by setting the gradient vector of Equation (85) with respect to w to zero, giving 
Y^E[e^(t)] = -2E[«(t)5(t)] + 2E[5(t)5"(t)]w, (87) 
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As per [27, 69, 73], the value of the weight veetor at time [t +1) is updated based on the 
method of steepest deseent as follows 

w{t + \) = w{t) + ^^[-Y^E[e^ (0]] , (88) 

where ju is the eonvergenee faetor. By substituting (t)] from Equation (87) into 

the above equation, we have 

= w(t) + (Off] 

= w(t) + jus{t^\jx{t^- (Off] 
whieh can be simplified to become: 

w(t + l) = w(t) + //5(t)e*(t), (89) 

where * represents complex conjugation and the instantaneous estimates given by 
[«(t)5(t)] and have been used instead of E'[«(t)5(t)] and 

^[^(0-^ (0] ’ respectively. As stated in [27, 69, 73], the term, £'[^ 5 (^) 5 ^ (t)] is 

referred to as the array correlation matrix Rss and implies a correlation between the 
signals received by the various array elements exists. The convergence factor // is related 

to the eigenvalues of ^(0^^ (0 • eigenvalues of ^(0^^ (0 widely spread, 
convergence will be slow. The value of // needs to be in the range shown in Equation 
(91) for convergence of the EMS algorithm [19, 73] 

0<//<^, (90) 

'^max 

where is the maximum eigenvalue of (t). 

The following is an iteration of the EMS adaptive beamforming algorithm. Eor 
our sensor grid, the central controller selects the number of array elements needed to 
form the desired beam. The central controller in effect creates an array, which is merely a 
subset of the total array. Eor example, a SIGINT grid that contains K elements randomly 
spread over an area, . The central controller will create a subset of A elements 
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randomly spread over an area, , where A<K and A 2 < A^. The eentral eontroller will 
then select F nodes contained within an area, A 2 , where PU A . 


[ /—xsinf^^)*cos( (f)+y sin(] 

^<{0,^) = ^ ^ ^ (91) 

where X is the wavelength of the signal in meters and x and y are the location of the 
individual antenna elements in meters. 0 is the elevation angle and (j) is the azimuth 
angle. 

The signals collected at these nodes are forwarded back to the central controller. 
The received signal is shown below. 

b{t) = A>{e4Ys{t) + n{t), (92) 

where 

s{t) = A[t)e \ (93) 

where / is the frequency of interest, t is the sample increment, is the sampling 
rate, y(^) is the phase and n{t^ is additive white Gaussian noise. Twice the frequency of 
interest was used as the sampling rate and five hundred data samples were used. The 
reference signal a{t^ needs to be closely correlated with ^(t). 


The central controller will then calculate the array weights w by summing each 
element over all possible values of 0 and (j). 


M N 


7.71 

,T, * y^f™sin(<?)cos(«»)+>'„„sin(6l)sin(. 


(94) 


where P = MxN the total number of nodes in the subset and is the complex weight 


» V^(^™sin(^t)cos(4)+>'„„sm((9jsin(4)) 

w =W e ^ 

mn mn 


(95) 


applied to the [m,nY element. The maximum value of 'F {0,^) occurs at 
[6,(P) = [6^,(Pq) , and the main lobe points towards ,^z^o) • 
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Each data sample of the signals collected from the ehosen nodes iteratively 
updates the corresponding array weights. The magnitude and angle of are 

stored. This eompletes an iteration of the adaptive beamforming solution. 

Should the eentral controller repeat this process using new data collected from the 
same P nodes or should the central eontroller request information from a different P 
node subset? Using the same P nodes could potentially deerease the network burden and 
reduee the number of required antenna elements. However, randomly selecting the nodes 
for eaeh subset eould distribute the power management aeross the network, as well as, 
lower the computational burden of the central controller. 

If the first option is ehosen then the eentral eontroller simply takes another set of 
data samples from each of the P nodes used in the previous calculation. The proeess 
described earlier is repeated and the final result is added to the previously stored result. 

If the seeond option is chosen then the central controller takes another set of data 
samples from a new subset of P nodes. Then the process deseribed earlier is repeated 
and the final result is added to the previously stored result. This magnitude is then plotted 
in dB against all possible values of azimuth and elevation. 

A third option is possible, as well. In this option the central controller will simply 
take another set of data samples from eaeh of the P nodes similar to the first ealeulation. 
However, in this ease the nodes will have moved in the intervening time interval. This 
will have the affect of creating a larger antenna field than truly exists. 

It will be shown that the second method is a mueh better way to calculate the 
weights. The first method is unable to overcome the grating lobes and provide a sufficient 
beam for direction of arrival determination. The second method beeause it is seleeting a 
different subset of nodes for eaeh iteration is able to overeome the grating lobes and 
provide a sufficient beam for direction of arrival determination with a minimal number of 
iterations. It will be shown that having mobile sensors will give an unparalleled 
surveillance capability with only a few sensors. However, it is assumed that these mobile 
sensors will be high eost and able to recharge their batteries at will. 
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The selection process is put into a general equation shown in Equation (96). How 
the central controller takes another set of data samples is what will be shown in 
subsequent paragraphs. Once the central controller selects more data the process 
described earlier is repeated and the final result is added to the previously stored result. 

( 96 ) 

(=1 


where L equals the total number of subsets created. Substituting Equation (95) into 
Equation (96) results in the following: 





m^=l n^=\ 


sin(^)cos{^i)+Tm.„. sin{6^)sin(^)j 


(97) 


where P = MxN the total number of nodes in each subset. This magnitude is then plotted 
in dB against all possible values of azimuth and elevation. This process is repeated until a 
desired solution has been found. 

1, Method One 

The selection using this method is described as follows. The central controller 
will randomly select a subset of P nodes from a total of A elements that are contained 
within an area, . As stated earlier, the overall array contains a total number of nodes. 


K, contained within and area, A^, where PU A< K and . These P nodes will 

form a subset of array elements that the central controller will use to calculate the array 
weights. The EMS algorithm adaptively produces weights for the formation of the beam. 
These weights are then used to calculate the magnitude and angle of the beam. Next the 
central controller obtains more data from the same subset of P nodes selected earlier. 
The new weights are then calculated and used to calculate the magnitude and angle of the 
beam in the desired direction. This result is then added to the previous calculations and 
then stored. This iterative process is continued for iterations where is used to 


represent Method One. Equation (97) is rewritten dropping the i dependency, to reflect 
the fact that the same P nodes are being used. This process is shown in Equation (98). 


= i, 


M N 

w e 

m=l n=l 


1.71 

j —sin(^)cos(^z^)+T^„ sm(^)sin(. 
A 


(98) 
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where P = MxN the total number of nodes in eaeh subset. This magnitude is then plotted 
in dB against all possible values of azimuth and elevation. This proeess is repeated until a 
desired solution has been found. MATLAB® eode developed in [73] was modified to 
eonduet the simulations in this seetion. This modified eode ean be found in Appendix I. 

MATLAB® eode in [73] was developed to randomly distribute sensors within a 
square grid. For ease of simulation within MATLAB the number of elements M was 
equal to the number of elements N. Onee this distribution of sensors was eompleted the 
eentral controller began to calculate the beam. As previously, stated the central controller 
continued to use the data received at the same P nodes. Several simulations were 
conducted to investigate the viability of this method. The results are as follows; 
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Figure 96. 25 Node array at a node density of 1 sensors/m^. 


The test parameters for this section assumed no errors and no node failures; 
elevation and azimuth 45°. Interferers were placed at 18° with elevation and azimuth at 
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45°. The number of nodes in the subset P. was seleeted to be 25. The number of 


iterations T, was varied between 10, 20 and 50 as shown in Table 18. The area, was 

varied between 25, 400, 2500 and 250000 m^. As was shown in Figures 32 - 34, as a 
given number of nodes are spread over an inereasing area the average distanee between 
sensors gets larger, therefore deereasing the node density. It ean be observed in Equations 
(24) - (27) as the frequeney goes up the beamwidth gets smaller and this was illustrated 
in Figures 31 (a) - (e). Therefore, by analyzing Equations (24) - (27) there are two 
faetors that have a big impaet on the beamwidth, one is the node density and the other is 
the frequeney of the signal-of-interest. These faetors will now be illustrated through 
several examples. The pertinent information for the simulations is shown in Tables 18 - 
20. The area, A 2 of 25 m is shown in Figure 96. 


Freq (MHz) 

Noise 

Power 

(pW) 

Area, 

A 2 

(m^) 

Sample 

Size 

#of 

iterations 

Li 

Number 

of 

Sensors 

per 

subset, 

P 

Node 

Density 

(sensors/m^) 

157/800/2400 

1 

25 

500 

10 

25 

1 

157/800/2400 

1 

25 

500 

20 

25 

1 

157/800/2400 

1 

25 

500 

50 

25 

1 


2 

Table 18. Adaptive beamforming utilizing Method One with a node density of 1 sensor/m . 
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Figure 97. The Array Pattern determined utilizing Method One with T/ =10 iterations for 
frequeneies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 


As ean be seen in the Figure 97 the beamwidth gets smaller as the frequeney goes 
up, whieh is eonsistent with the results shown in Figures 31 (a) - (e). This ean also be 
verified by observing the relationship of wavelength and node density to the beamwidth 
illustrated in Equations (24) - (27). Therefore for a given node density the elements are 
spaced at more wavelengths apart for the high frequencies. Also present at the higher 
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frequencies are a lot of side lobes that would need to be resolved if this node density was 
used. Ten iterations does not appear to be sufficient for this node density. 
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Figure 98. Array Pattern determined utilizing Method One with T; = 20 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 


As can be seen from Figure 98, 20 iterations provided a somewhat cleaner 
solution, but still insufficient to provide direction of arrival determination capability. 
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Figure 99. Array Pattern determined utilizing Method One with L/ = 50 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 

As can be seen from the figures, the gain of the mainbeam went from 25 dB for 
Li=10 iterations to 35 dB for Z;=50 iterations. Unfortunately, it can be seen that this 
method is not capable of mitigating the grating lobes. As the iterations are increased the 
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sidelobes have continued to increase and the mainlobe has continued to decrease. 
Therefore this method would not be capable of performing direction finding. 
Now the area, A 2 will be increased to 400 m as shown in Figure 100. 



2 

Figure 100. 25 Node array at a node density of 1.25 sensors/m . 


Freq 

(MHz) 

Noise 

Power 

(pW) 

Area, A 2 
(m') 

Sample 

Size 

#of 

iterations 

Li 

Number 
of Sensors 
per 

subset, P 

Node 

Density 

(sensors/ 

m^) 

157/800/24 

00 

1 

400 

500 

10 

25 

0.0625 

157/800/24 

00 

1 

400 

500 

20 

25 

0.0625 

157/800/24 

00 

1 

400 

500 

50 

25 

0.0625 


Table 19. Adaptive beamforming utilizing Method One with a node density of 0.0625 

sensor/m^. 
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Figure 101. Array Pattern determined utilizing Method One with L/ =10 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 


138 





Once again the side lobes are present at 800 MHz and 2400 MHz. Also it is worth 
noting that the solution has larger sidelobes than in the previous figures. This intuitively 
makes sense because with more distanee the effeet of the grating lobes beeomes more 
prominent. 
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Figure 102. Array Pattern determined utilizing Method One with T/ = 20 iterations for 
frequeneies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 
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Increasing the number of iterations has not reduced the side lobes at 800 MHz and 
2400 MHz. The mainbeam has increased to 25 dB, but the beamwidth is getting wider as 
the side lobes increase. 
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Figure 103. Array Pattern determined utilizing Method One with Z; = 50 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 

The mainbeam is now at 35 dB and the sidelobes have increased, which is similar 
to the case for the 25 square meter grid. Also the mainbeam beamwidth has continued to 
get wider. 
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Next the area, A 2 , was inereased to 2500 vcv as shown in Figure 104. 



2 

Figure 104. 25 Node array at a node density of 0.01 sensors/m . 


Freq 

(MHz) 

Noise 

Power 

(pW) 

Area, A 2 
(m^) 

Sample 

Size 

#of 

iterations 

Li 

Number of 
Sensors per 
subset, P 

Node 

Density 

(sensors/m^) 

157 

1 

2500 

500 

10 

25 

0.01 

800 

1 

2500 

500 

20 

25 

0.01 

2400 

1 

2500 

500 

50 

25 

0.01 


Table 20. Adaptive beamforming utilizing Method One with a node density of 0.01 

sensor/m^. 
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Figure 105. 


Array Pattern determined utilizing Method One with Z; =10 iterations for 
frequencies 157 MFlz (red), 800 MHz (blue) & 2400 MHz (green). 
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Figure 106. Array Pattern determined utilizing Method One with L; = 20 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 
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Figure 107. Array Pattern determined utilizing Method One with Z; = 50 iterations for 
frequencies 157 MHz (red), 800 MHz (blue) & 2400 MHz (green). 

The results for the 2500 square meter grid are analogous to the ones presented for 
the 400 square meter grid. Side lobes spikes are still present for 800 MHz and 2400 MHz 
at the 10 iterations solution. As the number of iterations increased the mainbeam went 
from 25 dB to 35 dB. The solution is noisier because the affect of the grating lobes has 

144 







become more of a problem as the distance between sensors increased. The sidelobes rise 
as the number of iterations increases. The beamwidth of the mainbeam is also getting 
wider. These last two facts are a direct result of the grating lobe problem. The algorithm 
is no longer selecting randomly from the entire sensor field therefore the build-up of the 
periodicities on the chosen elements becomes more prevalent as the number of iterations 
increases. 


Therefore it may be concluded that this approach is not adequate for direction 
finding particularly in a highly dense environment. In fact as more iterations are 
calculated the worse the direction finding performance becomes. Specifically, as the 
number of iterations increase the higher the side lobes become and the wider the 
mainbeam beamwidth becomes which means less directivity. In other words this 
methodology is not capable of mitigating the grating lobes. The amount of interference 
received from these side lobes would make it very difficult to localize a target of interest. 


In addition, the computational burden on the central controller and the energy 
burden on the individual sensors would be significant. For example, where total 
transmissions are defined as 

T^=hP, (99) 

where T, is the total number of iterations required for Method One to form the desired 

beam and P is the number of sensors involved in the process. Therefore, P is the 
number of transmissions per iteration in order to obtain a solution. If the cost per 
transmission is d joules / transmission , X\ joules is defined as the energy expended 
waiting to transmit. As stated before a node awaits for all other nodes to transmit prior to 
it transmitting. For convience, X\ is defined as a fraction of the energy to transmit the 

signal data. Xi joules is the energy expended to wake a node and assumes the worst case 


that all the nodes are in a sleep mode. For convience, Xi is also defined as a fraction of 
the energy to transmit the signal data. The total energy cost to form one beam in a given 
direction is defined as 


T^P=dT^+PX,+PX2 
= ^(TiP) + (P-l)(0.1^) + P(0.5^)’ 


( 100 ) 
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where F is in units of joules. For this method the energy eost is distributed over the same 
P subset of sensors resulting in an energy eost per sensor of 

S = ^F,+(0.1^)-^ + (0.5^), (101) 

where S is measured in units of joules/sensor. 

Note two of the examples diseussed earlier where the data is shown in Table 21. 
Three different frequeneies were ehosen for eomparison, and the node densities were all 
the same. In all eases the subset size was P=25 nodes. 


Freq 

(MHz) 

# of 

Trans., 
Tt 

Area, 

A 2 

(m^) 

Total 

Energy, 

Pv 

(joules) 

# of 
Iter,, 
Li 

Energy 

per 

sensor, 

(j/sensor) 

157 

/800 

/2400 

500 

25 

514.95 

20 

20.5965 

157 

/800 

/2400 

1250 

25 

1264.96 

50 

50.5966 


Table 21. Total number of transmissions and energy expended utilizing Method One. 


As ean be seen from Equation (99) that as F, inereases so does the total number 

of transmissions, therefore inereasing the total energy use per sensor as seen from 
Equation (101). Eor example, an inerease of F; from 20 to 50 inereased the total 
transmissions Fy/from 500 transmissions to 1250 transmissions. That is an inerease of 
total energy eost of 750(5 joules. To make matters worse inereasing the number of 
transmissions did not mitigate the grating lobes. Therefore this method is not suitable for 
beamforming and is espeeially not suited for direetion finding. 

2, Method Two 

Now it will be shown that randomly seleeting a different P subset of nodes per 
iteration from a large SIGINT array produees mueh better results. Using Equation (96) 
the eentral eontroller takes another set of data samples from a new subset of P nodes. 
Then the proeess deseribed earlier is repeated for L 2 iterations where L 2 is used to 
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represent Method Two. The final result is added to the previously stored result. This 
magnitude is then plotted in dB against all possible values of azimuth and elevation. This 
proeess is repeated until a desired solution has been found. 

Again the nodes were randomly distributed within a square grid. The number of 
M elements equals the number of N elements to faeilitate easier simulations. The software 
taken from [72] was modified to ereate the effeet of the eentral eontroller seleeting data 
from a different P subset of nodes for eaeh iteration. This modified eode ean be found in 
Appendix I. Above, the node density was varied by having the eentral eontroller seleet 
sensors that were distributed over different areas, A 2 . Another way to ehange the node 
density is by having more sensors, A, plaeed within a given area, A 2 , and was 
investigated. For an array oontainingzi > 75 sensors the area, A 2 , was ehanged between 
25, 2500 and 250000 m to eharaeterize the effeets on the beamforming proeess. 


Freq (MHz) 

Noise 

Power 

(PW) 

Area, A 2 
(m^) 

Sample 

Size 

#of 

iterations 

L 2 

Number of 
Sensors per 
subset, P 

Node 

Density 

(sensors/m^) 

157/800/2400 

1 

25 

500 

3 

>75 

>3 

157/800/2400 

1 

2500 

500 

3 

>75 

>0.03 

157/800/2400 

1 

250000 

500 

3 

>75 

>0.0003 


Table 22. Adaptive beamforming utilizing Method Two with varying node seleetion. 
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75 randomly chosen sensors from a node array with a node density > 3 

sensors/m^. 
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Figure 109. Array Pattern for a node density of >3 sensors/m for frequeneies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 

As can be seen in the plots the beamwidth gets smaller as the frequency goes up. 
Also present at the higher frequencies are a lot of sidelobe spikes that make the polar plot 
unreadable. This node density does not appear to provide satisfactory results. 
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Figure 110. 75 randomly chosen sensors from a node array with a node density > 0.03 

sensors/m^. 
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Figure 111. Array Pattern for a node density of >0.03 sensors/m for frequencies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 

The sidelobes in the Figure 111 (a) make it unreadable. Upon closer observation, 
it is noticed that the sidelobes are much higher than in Figure 109. This can be attributed 
to the fact that the distance has increased between the nodes, making the grating lobes 
more of a problem. It is also worth noting that the beamwidth of the mainbeam narrowed 
considerably as expected since the node density has decreased. 
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Figure 112. 75 randomly chosen sensors from a node array with a node density > 0.0003 

sensors/m^. 
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Figure 113. Array Pattern for a node density of >0.0003 sensors/m for frequeneies 157 

MHz (red), 800 MHz (blue) & 2400 MHz (green). 

The figure on the left still has too many side lobes present to be of any use. 
However, the figure on the right shows a very narrow beam. It ean be observed that the 
beamwidth of the mainbeam has deereased dramatieally. These results lead to the 
eonelusion that to narrow the beam, the further apart the sensors are spaeed the more 
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narrow the beam. A more in-depth diseussion on beamwidth will be presented in 
subseetion 4 of this chapter. Now an array size ofA> 125 sensors will be investigated. 
Again the area, A 2 , was altered between 25, 2500 and 2500000 m . In the next section 
since the overall number of sensors will be increased, this results in an increase in node 
densities. 


Freq (MHz) 

Noise 

Power 

(PW) 

Area, A 2 
(m') 

Sample 

Size 

#of 

iterations 

L 2 

Number of 
Sensors per 
subset, P 

Node 

Density 

(sensors/m^) 

157/800/2400 

1 

25 

500 

5 

>125 

>5 

157/800/2400 

1 

2500 

500 

5 

>125 

>0.05 

157/800/2400 

1 

250000 

500 

5 

>125 

>0.0005 


Table 23. Adaptive beamforming utilizing Method Two with a varying node. 
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Figure 114. 125 randomly chosen sensors from a node array with a node density > 5 

sensors/m^. 
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Figure 115. Array Pattern for a node density of >5 sensors/m for frequeneies 157 MFlz 

(red), 800 MFlz (blue) & 2400 MHz (green). 

These figures show a niee beamwidth formation for all three frequencies. With a 
few rogue sidelobes for 800 MHz and 2400 MHz frequencies, but much less than in the 
previous figures. 
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Figure 116. 
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125 randomly chosen sensors from a node array with a node density > 0.05 

sensors/m^. 
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Figure 117. Array Pattern for a node density of >0.05 sensors/m for frequencies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 


The results are as expected. The mainbeam beamwidth has been drastically 
reduced. The sidelobes appear much noisier. 
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Figure 118. 125 randomly chosen sensors from a node array with a node density > 0.0005 
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sensors/m . 
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Figure 119. Array Pattern for a node density of >0.0005 sensors/m for frequeneies 157 

MHz (red), 800 MHz (blue) & 2400 MHz (green). 


As can be observed from the figures, the mainbeam beamwidth has been reduced 
significantly. There is also an increase in sidelobes as the node density has been reduced. 
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Now an array size o^A>175 sensors will be investigated. The area, A 2 , is altered 

2 

between 25, 2500 and 2500000 m . This will have the affect of increasing the node 
densities at the respective areas. 


Freq (MHz) 

Noise 

Power 

(PW) 

Area, A 2 
(m') 

Sample 

Size 

#of 

iterations 

L 2 

Number of 
Sensors per 
subset, P 

Node 

Density 

(sensors/m^) 

157/800/2400 

1 

25 

500 

7 

>175 

>7 

157/800/2400 

1 

2500 

500 

7 

>175 

>0.07 

157/800/2400 

1 

250000 

500 

7 

>175 

>0.0007 


Table 24. Adaptive beamforming utilizing Method Two with a varying node. 
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Figure 120. 175 randomly chosen sensors from a node array with a node density > 7 

sensors/m^. 
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Figure 121. Array Pattern for a node density of >7 sensors/m for frequeneies 157 MFlz 

(red), 800 MFIz (blue) & 2400 MFIz (green). 


This figure shows a nice well defined beam with few rogue side lobes. As 
expected the mainbeam beamwidth decreases as the frequency goes up. 
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Figure 122. 175 randomly chosen sensors from a node array with a node density > 0.07 

sensors/m^. 
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Figure 123. Array Pattern for a node density of >0.07 sensors/m for frequencies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 
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Figure 124. 175 randomly chosen sensors from a node array with a node density > 0.0007 

2 

sensors/m . 


164 








Adfl()hve n«Rnilo(ming • FtaRrnMirltti ConpMisan 


157 MU/ 
800 MU/ 
2 4 GHj 


90 


20 


160 



0 


270 


iMe mntoi giid 


2 

Figure 125. Array Pattern for a node density of >0.0007 sensors/m for frequeneies 157 

MHz (red), 800 MHz (blue) & 2400 MHz (green). 
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As expected the mainbeam beamwidth for all frequencies has continued to narrow. 

Now an array size ofA> 250 sensors will be investigated. The area, A 2 , is altered 
between 25, 2500 and 2500000 m . This will have the affect of increasing the node 
densities at the respective areas. 


Freq (MHz) 

Noise 

Power 

(PW) 

Area, A 2 
(m') 

Sample 

Size 

#of 

iterations 

L 2 

Number of 
Sensors per 
subset, P 

Node 

Density 

(sensors/m^) 

157/800/2400 

1 

25 

500 

10 

>250 

>10 

157/800/2400 

1 

2500 

500 

10 

>250 

>0.1 

157/800/2400 

1 

250000 

500 

10 

>250 

>0.001 


Table 25. Adaptive beamforming utilizing Method Two with a varying node. 
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Figure 126. 250 randomly chosen sensors from a node array with a node density >10 

sensors/m^. 
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Figure 127. Array Pattern for a node density of >10 sensors/m for frequeneies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 
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Figure 128. 250 randomly chosen sensors from a node array with a node density > 1 

sensors/m^. 
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Array Pattern for a node density of >1 sensors/m for frequencies 157 MHz 
(red), 800 MHz (blue) & 2400 MHz (green). 








As expected the mainbeam beamwidth has been narrowed considerably. 
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Figure 130. 250 randomly chosen sensors from a node array with a node density > 0.001 

sensors/m^. 
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Figure 131. Pattern for a node density of >0.001 sensors/m for frequencies 157 MHz 

(red), 800 MHz (blue) & 2400 MHz (green). 


171 






It has been shown that randomly selecting a different Pi subset of nodes per 
iteration from a large wireless antenna array with a suitable node density produces much 
better results. Now it will be shown how this method also decreases the computational 
burden on the central controller and the energy burden on the entire sensor grid. This 
method assumes that elements contained within the area, A 2 , are only selected once, 
which is a reasonable assumption for a large array with equal probability of selection. 

The total transmissions required becomes 

7;=X/’=A, (102) 


where the subscript i in Pi is used to indicate that a different subset of nodes are chosen 
for each iteration. Again L 2 is use to represent Method Two. After modifying equation 
(100), for this methodology, which again assumes the worst case, that all nodes are a 
sleep and each node has to wait until the other sensors have transmitted, the total energy 
cost to form one beam in a given direction would be 
Tgp=dT^+Xx{^-^) + X2^ 


= ^£(/’) + (A-1)(0.1^) + A(0.5^) 
/=1 


(103) 


For this method, it is assumed that each node is only selected one time. This is a 
reasonable assumption if A is large and each node is equally likely to be selected. The 
energy cost is then distributed over all A sensors whereas in the first method the energy 
cost was spread over P sensors. This results in an energy cost per sensor of 

S=^ + (0.W)-5^+(0.55) 

£/; ■ (104) 

i=l 

S = ^ + (0.1^)-^ + (0.5^) 


This is roughly a factor T; reduction in the energy cost per sensor compared to the 
first method. The savings is even more significant when realizing that T; >> L 2 . 


The data for this experiment is shown in Table 26. The same three frequencies 
were chosen for comparison. For comparison only the energy expended to actually 
transmit the signal is used. Two different subset areas were chosen to demonstrate the 
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effects on the beamwidth of the main lobe. In all cases Pi = 25 nodes and L 2 = 10 
resulting in a total transmission cost of 250 transmissions. That is a reduction in total 
energy of 865(5 joules per beam. Not only has the overall energy savings increased, but 
the energy savings per sensor are quite significant. For the previous example the energy 
cost per sensor using the first method would be 50.596(5 joules/sensor, whereas for the 
second method the energy cost per sensor is (5 joules/transmission, which is a significant 
energy savings. 


Freq 

(MHz) 

# of 
Trans,, 
Tr 

Spacing, 
A 2 (m^) 

Total 

Energy, 

Pv 

(Joules) 

# of 
Iter., 
L 2 

Energy 

per 

sensor. 

157 

/800 

/2400 

250 

25 

399.93 

10 

\.65 

157 

/800 

/2400 

250 

2500 

399.93 

10 

1.63 


Table 26. Total number of transmissions and energy expended utilizing Method Two. 


As can be seen in the following figures, the mainbeam beamwidth has continued 
to get narrower as the central controller chose nodes that are further away from its own 
location. As expected as the main beam gain became narrower the side lobes increased. 
However, the main beam is large enough to mitigate the increase in side lobe level. 

The simulations clearly support what can be analyzed in Equations (24) - (27), 
that beamforming process is highly dependent on the node density and the frequency of 
operations. 

The results are as expected. Since the magnitude, after each iteration is added to 
the previous result and each iteration is allowing a new P number subset perspective the 
overall result is much better. The initial method while requiring less real sensors would 
be unable to perform the necessary SIGINT functions. 
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3, Beam width 

When conducting direction finding it is important to understand how to control 
the beamwidth of the beam. When in a search mode a wider beam would be preferred. 
After the signal has been localized it would be preferable to narrow the beam in order to 
enhance collection and/or facilitate offensive information operations. The approximations 
for beamwidth given earlier in Equations (24) - (27) are repeated here for convenience 
and shown as Equations (105) - (108) 


0 


elevation 


cos' (6»,) cos' (4) + sin' {(/),) 


0 


azimuth 


+ cos'(^/io)_ 


(105) 


(106) 


6^ = cos ' 


6 .. =cos ' 


cos6’„ -0.443 


cos6’„ +0.443 




[Length+ d) 

A 


^[Length+ d) 
where d=Length/N and N is the number of nodes. [3] 


(107) 


4 = COS^ 




cos -0.443 


cos + 0.443 




[Length+ d) 

A 

[Length+ d) 


(108) 


The solid angle and the directivity from Equation (28) and (29) are repeated for 
convenience and become Equations (109) and (110). 


^.4 ^elevation^azimuth ’ 


(109) 


on 


32400 
(degrees) ’ 


( 110 ) 
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As stated earlier an inverse relationship exists between the node density and the 
beamwidth. Likewise an inverse relationship ean be observed between beamwidth and 
direetivity. Speeifically, given a network with K number of nodes dispersed over an area, 
Ai, the node density ean be increased by selecting nodes that are dispersed within an area, 
A 2 where Aj <Aj to create the beam. Based on the size of the area, Aj, and the number of 
nodes, the beamwidth of the beam can be tailored to specific applications, specifically a 
wide beam for detection and a narrow beam for localization and/or collection. 

Figure 132 illustrates that when the number of sensors is increased, but the length 
of the array is held constant, i.e., the node density increases, the beamwidth increases 
resulting in a loss of gain, which can be seen in Equations (105) and (108). In general the 
more sensors that are added the directivity, D, i.e., gain, G, in a particular direction 
increases, but it is assumed that as more sensors are added the area is increased due to the 
increase in the physical area of the antenna. 


Beamwidth vs. Number of nodes, 157 MHz 



Figure 132. Beamwidth vs. Node Density where the blue asterisk represents measured 

values. Frequency is 157 MFlz. 
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Figure 132 demonstrates and is eonsistent with Equations (105) - (108) that a 
limit is reaehed on how mueh eontrol the eentral eontroller will have based on node 
densities and array distribution. 

Nonetheless, this is an interesting diseovery showing how the eentral eontroller 
ean tailor the beamwidths to speeifie funetions, e.g., seareh, loealization, eolleetion, 
offensive information operations, ete. When looking elosely at the node densities 
represented in the above Figure 132 it beeomes obvious that a limit on node densities has 
to be eonsidered. Nodes spaeed too elose together would potentially ereate a mutual 
eoupling problem that would need to be mitigated. More importantly, sinee it is 
envisioned that these sensors will need to be ineonspieuous and randomly deployed, too 
many sensors in a given area eould lead to a eompromise of intent. 

Now for a speeifie example, the software from Appendix 1 was modified to show 
the impact of sensor spacing and scanning a 360 degree arc. This modified software can 
be found in Appendix M. As seen in Tables 27 - 28 below a randomly distributed array 
containing K = 1000 sensors spread over an area of ^4, = 100000 formed the basis of 
our sensor grid. The frequency was kept constant at 157 MHz. 

Both examples will utilize the technique discussed in Chapter VI, Section B, 
subsection 2 to adaptively form the beams. The first example used sensors that were 
space within an area of = 400 that contained the central controller. A randomly 
chosen subset of P = 25 nodes was created for each iteration, to form the beam. For 
this example Fj = 10 • Focusing on the 60° and 45° beams from Figure 133, it can be 

observed that a 360° scan could be conducted at 15° increments with sufficient overlap to 
not leave gaps in coverage, but provide enough spread to limit the total number of beams 
that needed to be formed. Thus, this 360° scan can be accomplished with 24 beams. 
Increasing the increments would decrease the number of beams to be created but also 
leaves the potential for creating gaps in coverage. Decreasing the gaps provides more 
overlap, but it also requires creating more beams. 
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Elevation 

Angle 

0 

(radians) 

Azimuth 

Angle 

9 

(radians) 

Number 

of 

sensors 

per 

subset, P 

Area, 

A 2 (m^) 

71/4 

71/3 

25 

400 

71/4 

71/4 

25 

400 

7I/4 

71/6 

25 

400 

71/4 

Till 

25 

400 

7X/4 

n/S 

25 

400 

7X/4 

71/9 

25 

400 

7X/4 

Tl/lO 

25 

400 


2 

Table 27. Adaptive Beamseanning utilizing Method Two where, A 2 = 400 m 


Elevation 

Angle 

0 

(radians) 

Azimuth 

Angle 

9 

(radians) 

Number 

of 

sensors 

per 

subset, P 

Area, A 2 
(m^) 

7X/4 

71/3 

25 

40000 

71/4 

71/4 

25 

40000 

71/4 

71/6 

25 

40000 

71/4 

71/7 

25 

40000 

71/4 

71/8 

25 

40000 

7X/4 

71/9 

25 

40000 

71/4 

7c/10 

25 

40000 


2 

Table 28. Adaptive Beamseanning utilizing Method Two where, A 2 = 40000 m 


The second example used sensors that were spaced within an area of 
A^ = 40000 that encompassed the central controller. Now investigating the use of 
sensors that are located further away from the central controller it can be seen that the 
beam has narrowed considerably. In Figure 134 below it can be observed that a 
reasonable scan increment would be 2°/beam in order to provide a sufficient overlap. 
This would require 180 separate solutions to be computed. 
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Figure 133. Array pattern for a beamsean at 400 m with > 250 nodes utilizing Method 

Two. 
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Figure 134. Array pattern for a beamsean at 40000 m with > 250 nodes utilizing Method 

Two. 


The above results are consistent with the previous analysis that the frequency of 
the SOI, and node density have the greatest effect on beamwidth. The central controller 
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can vary the node density by seleeting larger areas to ereate the subsets. The subset 
ereation will be used to transition the beam from a signal seareh mode to signal 
proseeution mode. 

This seetion has demonstrated how to eontrol the beamwidth of the mainbeam. 
Also, shown was how a limit exists on how wide a beam ean beeome given a eertain node 
density. This limit far exeeeds the praetieal limit of the size of the sensor. Therefore, this 
researeh will use the array pattern displayed in Figure 133 for seareh. The array pattern 
displayed in Figure 134 will be used for a more foeused sustained eolleetion and/or 
offensive information operations posture. It is envisioned that both sensor grids used to 
ealeulate the beams will be subsets of a larger grid being eontrolled by the same beam 
eontroller. 

This ehapter has provided the reader an insight into time differenee of arrival and 
beamforming. The next ehapter will provide a more indepth analysis into the speed and 
energy tradeoffs between time differenee of arrival and beamforming 
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VII. DIRECTION FINDING METHODOLOGY 


The proposed direetion finding methodology for a SIGINT/IW sensor grid using a 
distributed antenna array is deseribed as follows. The SIGINT/IW wireless sensor 
network will eonsist of a large number of small sensor nodes that are densely deployed 
over an area of interest to aequire information about signals of interest. These sensor 
nodes eollaborate among themselves to form an ad-hoe network and disseminate the 
eolleeted target information to a eentralized gateway node, i.e. UAV, UUV, ete. The 
objeetive is to give national deeision makers or eombatant eommanders the eapability to 
have a relatively inexpensive solution for sustained SIGINT/IW or a rapidly deployable 
SIGINT/IW alternative. 

A two-tiered hierarehieal elustering sensor network arehiteeture is assumed. 
Distributed array establishment and eoordination will be aehieved by using elustering 
algorithms to form a hierarehieal arehiteeture. In this researeh it is assumed that the 
elustering funetions to ereate a sensor grid have been eompleted. It is assumed that 
perfeet node loeation has been determined and that this data is available to the eentral 
eontroller. A eentral eontroller (i.e. UAV, UUV, ete.) is then tasked with maintaining 
frequeney, phase and data synehronization among the remaining nodes (or seeondary 
nodes) within the eluster. The primary node aets as a referenee antenna by gaining initial 
intereept and frequeney-of-interest (FOI) determination. The primary node will then 
eoordinate the DOA determination. Onee the DOA of the SOI has been determined the 
primary node will eoordinate the beamforming proeess to direet an antenna beam towards 
the SOI in an adaptive manner in order to maximize eolleetion or faeilitate offensive 
information operations. 

Now a more detailed diseussion will be offered into how best to eonduet direetion 
finding with the unique arehiteeture of the wireless SIGINT/IW antenna network. This 
researehed foeused on two eompeting methodologies for direetion-fmding within a 
wireless antenna array network with the speeifie foeus being on rapid resolution of a line- 
of-bearing. The first methodology is adaptive beamforming that performs a beamsean 
similar to a phased-array radar. The seeond methodology uses a eombination of TDOA 
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and adaptive beamforming. The TDOA lines-of-bearing are calculated using a 
combination of the Cross-Ambiguity Function to calculate the actual time difference and 
the Newton-Raphson technique to compute a geolocation. In both cases the adaptive 
beams are formed using a LMS algorithm. In this chapter, these methods have been 
adapted to the unique random distribution of the wireless sensors. This comparison is 
done to find the most expeditious way to calculate a solution and determine the angle-of- 
arrival. 

A. ENHANCED COLLECTION METHODOLOGY (ECM)-l 

Now a look at the time required to perform a combination TDOA and 
beamscanning solution. The same set-up requirements from Section A are assumed. 

Step 1. The central controller acts as a reference antenna by gaining initial 
intercept and FOI determination. 

Step 2. Based on the analysis in Chapter VI, Section A, subsection 2 the central 
controller will then determine the TDOA collected from the fl sensors furthest from its 
location. Each sensor forms a two-node collector pair with the central controller. The 
central controller will utilize the cross-ambiguity function presented earlier to calculate 
the TDOAs between each of these fl sensors and the central controller. 

Step 3. The central controller will then use the modified Newton-Raphson 
technique to localize the signal-of-interest to a specified quadrant. 

Step 4. The central controller will then coordinate the DOA determination using 
beamscanning. The central controller will create a subset of A elements randomly spread 
over an area, A 2 , where A <K and A 2 <Ai. The central controller will then select Pi nodes 
contained within an area, A 2 , where P « A. Following the algorithm in Chapter VI, 
Section B, subsection 2, the central controller will create L subsets utilizing a different Pj 
nodes for each iteration to create the appropriate beam for a given {6^,(I)q). 

Step 5. As stated above, in Section A, based on the results of this search then the 
beam controller can randomly choose sensors that are dispersed within a larger area to 
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further narrow the beam into the target direetion in order to maximize eolleetion or 
faeilitate offensive information operations. 

Now the total time required to do eombination TDOA and beamforming ean be 
defined as follows: 

^TDOAI'V ~~ ^TDOA ^BF ’ (HI) 

where is the time to required to loealize the target and is the time required to 

form one beam. After substituting Equation (78) into Equation (111) and eombining like 
terms the total time required ean be rewritten as 

^TDOAI'V ~ + Os' + Otj + 0^2 + tip + Atj + At2 5 (112) 


where h is the time required to transmit and is the time required to bring a sensor out 

of sleep mode. Equation (112) assumes the worst ease that all the nodes need to be 
brought out of a sleep mode. In addition, it assumes that the nodes used in the TDOA 
ealeulations are not used in the beamforming ealeulations. 

The total energy to enhanee eolleetion using a eombination of TDOA and 
beamforming is defined as 


E =r 

^ TDOAI'V ^ TDOA 


(113) 


where again is defined in Equation (76) as the energy required for TDOA and E^^^^ 

is defined in Equation (103) as the energy required to form one beam using the multiple 
random seleetions method. After substitution, Equation (113) ean be rewritten as follows: 

^TDOAI'V ~ Q^+(Q-i)Y;ri+H;r2+^Z^+(A-i);ri+A;r2 

i=\ 

L 2 ’ 

=Q^ + (Q-1)Y(0.1^) + Q(0.5^) + ^2/)+(A-1)(0.1^) + A(0.5^) 

/=! 

(114) 


B, ENHANCED COLLECTION METHODOLOGY (ECM)-2 

The proposed approaeh for a wireless antenna sensor grid using beamseanning 
only is deseribed as follows. Eor our sensor grid, the eentral eontroller seleets the number 
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of array elements needed to form the desired beam. The central controller in effect 
creates an array, which is merely a subset of the total array. For example, a SIGINT grid 
that contains K elements randomly spread over an area, Aj. 

Step 1. The central controller acts as a reference antenna by gaining initial 
intercept and FOI determination. 

Step 2. The central controller will then coordinate the DOA determination using 
beamscanning. The central controller will create a subset of A elements randomly spread 
over an area, A 2 , where A <K and A 2 <Ai. The central controller will then select Pi nodes 
contained within an area, A 2 , where P « A. Following the algorithm in Chapter VI, 
Section B, subsection 2, the central controller will create L subsets utilizing a different Pi 
nodes for each iteration to create the appropriate beam for a given {6^,(/)q). 

Step 3. The beam controller will repeat Step 2 until a 360 degree scan has been 
completed. 

Step 4. Once the SOI has been localized to a specific quadrant then the central 
controller will repeat step 2 , except the central controller will use nodes that are dispersed 
over a larger area, ^4^ to create the Pi node subset. This will allow the central controller to 
create a narrower beam to localize the SOI to a specific sector of the quadrant. More 
sensors can be incorporated into the adaptive beamforming solution if higher gain is 
required. 

Step 5. Based on the results of this search the beam controller can randomly 
incorporate sensors that are dispersed within an even larger area, A4, to further narrow the 
beam into the target direction in order to maximize collection or facilitate offensive 
information operations. 

Now the total time required to do beamscanning only can be defined as follows: 

hs = ^4 + + Atj, (115) 

where ttp is the time to required to form one beam, as defined earlier, and 7 is the number 
of beams required to conduct a 360 scan. ?/ is a function of the solid angle of the beam 
and the user defined overlap desired, h is the time to transmit and is the time required 
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to bring the node out of a sleep eyele. Equation (115) assumes the worst ease therefore 
providing an upper bound on the time required to perform ECM-2. Equation (115) 
assumes that a transmitting node will have to wait until all the other nodes are finished 
transmitting. Eikewise Equation (115) assumes that all nodes partieipating in ECM-2 will 
only need to be brought out of a sleep mode one time. 

The total energy to enhanee eolleetion using beamseanning only is defined as 

^bs=ri{5T^) + ri{h-\)x,+^X2 
= ;7(^7;) + ;7(A-1)(0.1^) + A(0.5^)’ 


where rj is the number of beams that need to be formed. Equation (116) is similar to 
Equation (103) exeept that multiple beams need to be formed. After substitution, 
Equation (116) ean be rewritten as follows: 


r 


bs 



V -1 J 


+ ;7(A-1)(0.1^) + A(0.5^), 


(117) 


Equation (117) assumes the worst ease seenario in order to provide an upper 
bound on energy usage. Equation (117) assumes that energy is spent waiting to transmit 
and that eaeh node had to be brought out of a sleep eyele. The first term in Equation 
(117) aeeounts for the number of beams to be formed. The seeond term in Equation (117) 
aeeounts for the energy expended waiting to transmit. Einally, the third term in Equation 
(117) is the energy expended to bring the partieipating nodes out of a sleep mode. It is 
assumed that the nodes will only need to be brought out of a sleep mode one time. 

Now a look at the eombination of TDOA and beamforming. 


C. COMPARISON 

The energy eost using the eombination of TDOA and beamforming will always be 
less than the beamseanning only method as long as the following eriteria is met: 
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r <r 

^ TDOAI'V — bs 


Q^ + (Q-1)Y(0.1^) + Q(0.5^) + ^ +(A-1)(0.1^) + A(0.5^) 

V '=1 ) 

f L, \ 

<7lS ^P. +;7(A-1)(0.1^) + A(0.5^) ,(118) 

V i=i y 

^2 

QY + QY(0.1^)-Y(0.1^) + Q(0.5^)<(;7-1) J]/’ + (A-1)(;7-1)(0.1^) 

i=l 

r,^o^+QY(0.1^)-Y(0.1^) + Q(0.5^)<(;7-l)7;+(A-l)(;7-l)(0.1^) 


Upon close inspection of Equation (118), it can clearly be seen that ECM-1 will 
always use less energy than ECM-2 as long as the number of transmissions to perform 
the TDOA plus the additional factor associated with the energy spent waiting to transmit 
and the energy required to bring a sensor out of a sleep cycle is less than the number of 
transmissions required to conduct the 360 beamscan minus the target beam. The target 
beam will be formed in both cases. 

Additionally, the combination of TDOA and beamforming will always be faster 
as long as the following criteria is met; 


^TDOAI'V — hs 


O^Y + Pis + OYtj + 0^2 + tvp + Atj + Atj ^ + At2 

Q^Y + Q^ + QYtj +0^2 - (^“1)4 

^TDOA 


,(119) 


Equation (119) shows that ECM-1 will be faster than ECM-2 will be faster as 
long as the time required to form the beams, less the target beam, plus the amount of time 
spent waiting to transmit and the time required to bring a sensor out of a sleep mode is 
greater than the time to determine the target line-of-bearing the combination of TDOA 
and beamforming will always be faster than beamscanning only. 

Eocusing on the 60° and 45° beams from Eigure 133, it can be observed that a 
360° scan could be conducted at 15° increments with sufficient overlap to not leave gaps 
in coverage. This can be accomplished with ;; = 24 beams. Using the beamwidth from 
Eigure 134 it would require 6850 time units to scan 360°. 
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The energy and time cost are dependent on the number of beams that need to be 
formed for beamscanning only. As stated earlier, the number of beams required would be 
based on the solid angle of the beam and the desired overlap of the beams. Table 29 
depicts the energy cost and times for various ?/. The t ic-toc command in MATLAB® 
was used to determine that t<p = 25 time units were required to create a beam, e was 
determined to be equal to 0.1 time units and g was determined to be 2 time units per CAP 
estimate for each two-node collector pair. These times would vary based on hardware and 
software implementation. Additionally, it is expected that these times would be much 
lower in a real world implementation. From Table 26, the energy to form one beam is 
399.9(5 joules. From Equation (76), the total energy to do TDOA is 45(5 joules, where Q = 
4 sensors. 


# of beams, r\ 

Tbs, (joules) 

r tdoa/t? 

(joules) 

tbs, (s) 

tTnOA/T, (s) 

10 

28745 

444.95 

3000 

612.4 

24 

6722.65 

444.95 

6850 

612.4 

30 

83725 

444.95 

8500 

612.4 


Table 29. Energy and Time Cost comparison between ECM-1 vs. ECM-2. 


It can be seen from the results shown in Table 29 that the combination of TDOA 
and beamforming in most practical cases will always be faster and more energy efficient 
than the beamscanning only method. 

Once the correct beam has been formed then the determination on whether to 
continue prosecution will need to be made. In addition, using the knowledge gained on 
how to affect the beamwidth of the beam, a narrower beam can be created for sustained 
collection and/or offensive information operations. Figure 135, shows a schematic 
diagram of the information flow in the SIGINT/IW wirelessly distributed antenna system. 
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Figure 135. Schematic diagram of the information flow in a SIGINT/IW wirelessly 
distributed antenna system. *Note - Synchronization implies both timing and 
local oscillator synchronization. Initial Conditions: time, t = 0. The following 
assumptions apply: All participating sensors are on; All sensors are 
interconnected and synchronized (timing and local oscillator); All nodes have 
connectivity with the central controller node; Sensor Receivers monitor a 
specified band, i.e. AD8347 (800 MHz - 2.5 GHz). 


In this chapter a methodology for rapid resolution to determine lines of bearing in 
a wireless antenna array network has been presented. Two competing methodologies 
were discussed. One method used beamscanning exclusively and the other method used a 
combination of TDOA and beamscanning. The metric for determining the best method 
was the time required to perform each method. It was shown that the combination TDOA 
and beamscanning method is far superior to the beamscanning only method. 
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VIII. SUMMARY 


In this research, a novel approach for conducting signals intelligence from a 
distributed network of wireless nodes was developed. The primary objective of this 
research is enhancing collection in a specified target direction. The focus was on 
enhancing RF signal collection through beamforming in a wireless antenna system. A 
variety of techniques were investigated with a combination of TDOA and adaptive 
beamforming showing the most promising results. The combination of the two techniques 
dubbed Enhanced Collection Methodology (ECM) - 1 provides rapid acquisition, 
localization, collection and/or the option of conducting offensive information operations, 
against a particular SOI. 


A, CONCLUSIONS 

In the development, of ECM-1 two specific areas of conducting RE collection 
were addressed in this research. One is the time required to determine the target direction 
and the time required to form the beams. The other is the energy consumption involved in 
developing these solutions. Two competing enhanced collection methodologies (ECM) 
were developed and researched, ECM-1 and ECM-2. ECM-1 uses a combination of time 
difference of arrival (TDOA) and adaptive beamforming. ECM-2 uses adaptive 
beamforming that performs a beamscan similar to phased-array radars. 

Additionally, two competing methods for forming the beams were created. 
Method One and Method Two. Method One created a subset of sensors from a randomly 
distributed antenna array. This method uses data exclusively from these elements and 
increases the number of iterations until the desired beam is formed. Method Two creates 
a new subset of sensors, from a randomly distributed array, for each iteration, until the 
desired beam is formed. 

Analytic expressions were derived that proved that Method Two provides better 
power management across the sensor network. Additionally, Method Two was able to 
mitigate the grating lobes that are present in arrays where the antenna element spacing is 
greater that half-wavelength. 
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Another important discover within this research was the ability of the central 
controller to control the beamwidth, therefore the gain of the main beam. By randomly 
selecting nodes that were distributed over different areas the central controller, in effect 
manipulates the node density. Outside of the frequency of the tartget signal the node 
density was shown to have the greatest affect on the beamwidth of the main beam. The 
significance of this discovery is that no additional energy burden was placed on the 
network. 

Analytic expressions, for energy consumption and time to determine a line of 
bearing solution, were derived to compare ECM-1 to ECM-2. The energy consumption 
for ECM-1 is based on the energy required to determine the line of bearing to the target 
using time difference of arrival plus the energy required to form the beam. This method 
takes advantage of the spatial diversity of the nodes. Multiple estimates taken from a 
select number of nodes located furthest from the central controller reduces the number of 
nodes required to determine the target direction. By reducing the number of participating 
nodes the overall energy consumption of the network is reduced. 

The energy required for ECM-2 is based on the requirment to form multiple 
beams in order to provide sufficient coverage for a given area. Beamforming was shown 
to require a lot of energy. Thus having to form multiple beams placed a tremendous 
energy burden on the network and was shown to be too costly. 

The time required to conduct ECM-1 is based on the time required to determine 
the target line of bearing through time difference of arrival plus the time required to form 
the beam in the target direction. This method also takes advantage of the spatial diversity 
of the individual nodes. Multiple estimates taken from a select number of nodes located 
furthest from the central controller is all that was required to determine a line of bearing 
to the target. After the line of bearing to the target has been determined then only the one 
beam needs to be formed. 

In comparison, the time needed to conduct ECM-2 is based on the requirment to 
form multiple beams in order to search a given area. The number of beams required to 
search a given area is a user defined requirement that is based on the beamwidth of the 
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main beam. Beamforming was determined to be a time consuming evolution and the 
requirement to form multiple beams to determine the target direction proved too costly. 

ECM-1 was shown to be far superior to ECM-2 in both energy consumption and 
the time required to formulate a solution. 


B, AREAS OF FUTURE RESEARCH 

A few potential follow-on research areas will now be discussed. One area of 
research is the ability of this methodology to prosecute mobile targets. Another research 
area is the ability of this network to work against wideband signals. ECM-1 should be 
robust enough after some minor modifications to address these new requirements. 

Characterizing the types of traffic will be crucial. This dissertation assumed that 
sufficient buffer space was resident at the individual nodes. Understanding the type of 
traffic would be critical in selecting the appropriate medium access control scheme, as 
well as, buffer sizes at the individual nodes. If buffers overflow then critical data needed 
to perform the beamforming would be lost. Selecting the appropriate MAC will be 
critical in insuring that sensors can access the network and transmit their data reliably 
with a minimum of retransmissions before their buffers overflow. It is expected that the 
same data will be arriving at the nodes at slightly different time intervals, therefore a 
TDMA scheme would probably be the best option. 

Other issues are the synchronization of the nodes. It was assumed in this 
dissertation that all sensor nodes were synchronized. It is believed that the 
synchronization demands will be such that the nodes will need to have constant 
communications with the central controller. If this is indeed the case then a separate 
channel will probably be needed to accommodate this requirement. Both the data and 
synchronization channels will need to be at a minimum low probability of intercept (EPI) 
channels to prevent compromising the intent of this network and exposing it to various 
attacks. It would be optimal to have the channels be low probability of detection (EPD) 
signals, as well. 

Due to the constant transmissions expected to keep the sensors synchronized 
research into power consumption will be important. The synchronization methods 
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discussed in this research are power intensive and assume that the nodes are connected to 
a reliable energy source, i.e., a ship’s power system. Transmissions are one of the highest 
energy burdens on a sensor. Therefore, innovative solutions will be needed for the type of 
modulation scheme designed for synchronization. 

A more intense study of the sample size should be investigated. In this research, 
500 samples of data were used for calculating a solution. Reducing the sample size would 
reduce the transmission length and the loading of the network. 

The use of mobile nodes will also merit further study, as well. In this dissertation, 
all nodes were assumed to be static. Having static nodes required a large number of 
nodes. If mobile nodes are used then a smaller number of nodes may be required. 

However, what will be the impact on the beamforming process of the samples 
being taken at different instances in time? This will no doubt have a deleterious effect on 
the formation of the beams, but the question becomes how much. 

Initial investigations show promising results. The direction finding method is 
similar to the first method discussed in Chapter VI. The central controller will select a 
subset of mobile nodes, P that are contained within an area, . These P nodes will 

form a subset of array elements that the central controller will use to calculate the array 
weights. The LMS algorithm adaptively produces weights for the formation of the beam. 
These weights are then used to calculate the magnitude and angle of the beam. Next the 
central controller obtains more data from the same subset of P nodes selected earlier, 
except unlike in the first method these nodes have moved. Therefore, the new data will be 
taken from a different location. The new weights are then calculated and used to calculate 
the magnitude and angle of the beam in the desired direction. This result is then added to 
the previous calculations and then stored. This iterative process is continued for L 
iterations. Equation (94) is rewritten dropping the i dependence, to reflect the fact that the 
same P nodes are being used. This process is shown in Equation (98). The only 
difference will be in the beamforming equation. Equation (97) will be modified to 
account for the z component of our sensors as shown in Equation (120). 
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(120) 


[ J—*(xsm( ^^)*cos(^)+ V sin(^?)sin{^)+z cos(1 

W{0,0) = e^ 

where A is the wavelength of the signal in meters; x, y, and z are the loeation of the 
individual antenna elements in meters; and 0 is the elevation angle and (j) is the azimuth 
angle. This magnitude is then plotted in dB against all possible values of azimuth and 
elevation. This proeess is repeated until a desired solution has been found. MATLAB® 
eode developed in [73] was modified to eonduct a few simulations. This modified code 
can be found in Appendix K. 

Something very interesting was discovered. If the cos(^) term in the exponent 

was dropped the beamwidth of the main lobe can be increased, thereby increasing the 
detection region in a search mode. 

This method may show the viability of using a few nodes and being able to gain 
some of the same benefits that a larger array would provide. Certainly, these nodes would 
be more expensive than their static counterparts. This method could be applied to any 
type of mobile sensor, e.g. UAV, UUV, etc. 
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APPENDIX A. BRUTE FORCE 


% This is the brute_force program for beam force technique. 

"6 

% Can use either stored Array or create an array on the fly, position 
of 400 elements and M*N1=400; 

clear 

%load('Array'); 

M=20; %size(x_err); 

N1 = M; 

L = 500; 
error2= 0.0; 

[x,x_err,y,y_err,distance,midpointX,midpointY]=Random_LOC_DOA(M,N1,L,er 
ror2) ; 

%x=x_err; 

%y=y_err; 

figure(l), plot(x,y, ' * ' ), grid on, hold on 
N=M*N1; 

x_errl=reshape(x_err,N,1)'; 
y_errl=reshape(y_err,N,1)'; 

fc=2400e6;k=2*pi/(3e8/fc);rad=pi/180; % element 

frequency and constants 

r=sqrt(sum([x_err1.*x_err1;y_err1.*y_err1]));[min_r ref]=min(r); % 

select reference element closest to origin 

phase_ref=mod(k*2*(r-min_r)/rad,360); % generate 

reference phase 

phase=zeros(1,N); % initialize 

phase-shifter 
%amp=l./(r/min_r); 
amp=ones(1, N) ; 

time=l; %initialize 

time tag 
for 1=1:N 

element_start(i)=time; 

error(time)=amp(i)*exp(j*(phase(i)*rad))*exp(j*k*2*r(i))- 
exp(j*k*2*r(ref)); 

errorp(time)=mod((phase(i)*rad+k*2*r(i)-k*2*r(ref))/rad, 3 60); 
if errorp(time)>180, 

errorp(time)=errorp(time)-360; 

end 

while ~((abs(error(time))<=0.2)||(time-element_start(i))>14) 
phase(i)=phase (i)+22.5; 
time=time+l; 

error(time)=amp(i)*exp(j*(phase(i)*rad))*exp(j*k*2*r(i))- 
exp(j*k*2*r(ref)); 

errorp(time)=mod((phase(i)*rad+k*2*r(i)-k*2*r(ref))/rad, 360) ; 
if errorp(time)>180, 

errorp(time)=errorp(time)-360; 

end 

end 

time=time+l; 

end 

figure(2) 
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%generate plot 



hold on, grid on 
for i=l:N 

a=element_start(i); 
if i==N, 

b=length(error) ; 

else 

b=element_start(i+1)-1; 

end 

if rem(i,7)==1, 

plot((a:b),errorp(a:b),'b') 
elseif rem(i,7)==2, 

plot((a:b),errorp(a:b),'g') 
elseif rem(i,7)==3, 

plot((a:b),errorp(a:b),'r') 
elseif rem(i,7)==4, 

plot((a:b),errorp(a:b),'c') 
elseif rem(i,7)==5, 

plot((a:b),errorp(a:b),'m') 
elseif rem(i,7)==6, 

plot((a:b),errorp(a:b),'y') 

else 

plot((a:b),errorp(a:b),'k') 

end 

end 

xlabel('Number of iterations') 
ylabel('Phase error (deg)') 
figure(3) 

final_errorp=mod(phase_reffphase,360); 
for i=l:length(final_errorp) 
if final_errorp(i)>180, 

final_errorp(i)=final_errorp(i)-360 

end 

end 

plot(final_errorp,'xr') 
xlabel('Element number') 
ylabel('Final Phase error (deg)') 
grid on 

axis ( [0 500 -40 40]) 
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APPENDIX B. BEAM TAGGING 


% This is the beam_tagl program for beam tagging technique. 

"6 

% Used stored Array created for comparison with brute_force, position 
of 400 elements and M*N1=100; 

clear 

load('Array'); 

M=size(x_err); 

N=prod(M); 

x_errl=reshape(x_err,N,1) 
y_errl=reshape(y_err,N,1) 

fc=3e8;k=2*pi/(3e8/fc);rad=pi/180; 
frequency and constants 

r=sqrt(sum([x_err1.*x_err1;y_err1.*y_err1]));[min_r 
select reference element closest to origin 
phase_ref=mod(k*2*(r-min_r)/rad,360); 
reference phase 
phase=zeros(1,N); 
phase-shifter 

% Effects of amplitude with distance 
%amp=l./(r/min_r); 
amp=ones(1, N) ; 

time=l; %initialize 

time tag 
command=[2]; 
for i=l:N 

element_start(i)=time; 
command=[command 0]; 

while command(length(command))tcommand(length(command)-1)~=0, 

F_pos(time)=amp(i)*exp(j *((phase(i)+90)*rad+k*2*r(i)))+exp(j*k*2*r(ref) 

); 

F_neg(time)=amp(i)*exp(j*((phase(i)- 
90 ) *rad+k*2*r(i)))+exp(j*k*2*r(ref)); 

if abs(F_pos(time))<abs(F_neg(time)), 

phase(i)=phase(i)-22.5;command=[command -1] ; 

else 

phase(i)=phase(i)+22.5;command=[command 1] ; 

end 

errorp(time)=mod((phase(i)*rad+k*2*r(i)-k*2*r(ref))/rad,360); 
if errorp(time)>180, 

errorp(time)=errorp(time)-360; 

end 

time=time+l; 

end 

end 

figure(1) %generate plot 

hold 

for i=l:N 

a=element_start(i); 


% element 
ref]=min(r); % 

% generate 
% initialize 
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if i==N, 

b=length(errorp) ; 

else 

b=element_start(i+1)-1; 

end 

if rem(i,7)==1, 

piot((a:b),errorp(a:b),'b') 
elseif rem(i,7)==2, 

piot((a:b),errorp (a:b), 'g') 
elseif rem(i,7)==3, 

piot((a:b),errorp(a:b),'r') 
elseif rem(i,7)==4, 

piot((a:b),errorp (a:b), 'c') 
elseif rem(i,7)==5, 

plot((a:b),errorp(a:b) , 'm') 
elseif rem(i,7)==6, 

piot((a:b),errorp (a:b), 'y') 

else 

piot((a:b),errorp (a:b), 'k') 

end 

end 

xlabel('Number of iterations') 
ylabel('Phase error (deg)') 
figure(2) 

final_errorp=mod(phase_reffphase, 360); 
for 1=1:length(final_errorp) 
if final_errorp(1)>180, 

final_errorp(1)=final_errorp(1)-360; 

end 

end 

plot(final_errorp,'xr') 

xlabel('Element number') 

ylabel('Steady Phase error (deg)') 

grid 

axis ( [0 500 -90 90]) 

%figure (7) 

%for i=l:N 

% if i==N, 

% count_tag(i)=time-element_start(i); 

% else 

% count_tag(i)=element_start(i+1)-element_start(i); 

%end 

%end 

%plot(count_tag-l, 'xr' ) 

%xlabel('Element Number') 

%ylabel('Number of iterations required') 

%figure (8) 

%hist(count_tag-l) 

%xlabel('Number of iterations required') 

%ylabel('Number of elements') 
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APPENDIX C. CAF I 


function [TDOA, FDOA] = CAF(SI, S2, Max_f, fs, Max_t); 

9 - -k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

% CAF takes as inputs two sampled signal vectors (SI & S2) in analytic 
% signal format, the maximum expected FDOA in Hertz (Max_f), the 
% sampling frequency used to generate SI & S2 (fs), and the maximum 
% expected TDOA in seconds (Max_t). The function then utilizes 
% Stein's method in [1] to compute coarse estimations of TDOA and 
% FDOA between SI & S2. Finally, "fine mode" calcualtions are made 
% to compute the final TDOA and FDOA, which are returned to the 
% user via the output arguments. 

% Written by: LCDR Joe J. Johnson, USN 
% Last modified: 17 September 2001 
% Modified by CDR Mickey Batson, USN 
% Last modified: 01 October 2006 

2 - -k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

clc; 

N = length (SI); 

51 = reshape(SI,N,1); % Ensures signals are column vectors due to 

52 = reshape(S2,N,1); % Matlab's better efficiency on columns 
Sl_orig = SI; % Want to preserve original input signals 
S2_orig = S2; % for later use; SI & S2 will be 

% manipulated in the fine mode below. 

% The following while loop ensures that the sub-block size, Nl, is 
% large enough to ensure proper resolution. If Max_f/fs*Nl were 
% less than 1, then the Freq calculated at the end would always be 
% + or - 1/Nl! 2^19 = 524288 is about the limit for efficient 
% processing speed. 

Nl=1024; 

while (Max_f/fs*Nl < 2) & (Nl < 2^19) 

Nl = 2*N1; 
end 

N2=Nl/2; 

if Nl > N % For cases where resolution calls for 

51 = [SI;zeros(Nl-N,1)]; % a sub-block size larger than the 

52 = [S2;zeros(Nl-N,1)]; % signal vectors, pad the vectors with 
N = Nl; % zeros so that they have a total of 

end % Nl elements. 

% Want magnitude of Max_f, since +&- will be used below 
Max_f = abs(Max_f) ; 

Number_of_Blocks = length (SI)/Nl; % Number of sub-blocks to break 
% the signal into 

Min_v = floor(-Max_f/fs*Nl); % Smallest freq bin to search 
Max_v = -Min_v; % Largest freq bin to search 
v_values = Min_v : Max_v; % Vector of all bins to search 
Max_samples = Max_t * fs; % Maximum number of samples to search 
% Finds max number of block shifts (q) that must occur for each 
% R and v below. 

If Max_samples > N2 

q_max = min(ceil((Max_samples - N2)/Nl),Number_of_Blocks-l); 
else q_max = 0; 
end 
x=0; 
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divisors = Number_of_Blocks:-1:1; % Used to scale "temp" below... 

9 - -k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

% COARSE MODE computations. 

2 - -k-k-k-k-kk-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-kk-k-k-k-k-k-k-k-k-k-k-k-k-k-k-kk-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

for V = 1:length(v_values) 

temp (1: N1, 1: q_max+l) =0 ; % Initializing - saves time.... 
for R = 0:Number_of_Blocks-l 

% tempi is the EFT of the R'th block of SI, shifted by "v" bins. 

Tempi = fftshift(fft(SI(1+R*N1 : N1*(R+1)))); 
tempi = shiftud(tempi,v_values(v),0); 
for q = 0:q_max 

% R+q cannot exceed the number of sub-blocks 

if R + q > Number_of_Blocks-l break 

end 

% EFT of the (R+q)'th block of S2 

temp2 = f ftshift (f ft ( [S2 (1+(R+q) *N1 : N2 + Nl* (R+q) ) ;... 
zeros(N2,1)])); 

% Multiplies tempi & temp2, FFTs the product, then adds to 

% previous values for the same value of q (but different R) 

temp(:,q+l) = temp(:,q+l) + ... 

abs(fftshift(fft(tempi.*conj(temp2)))); 

end 

end 

% Each value of q was used a different # of times, so they must be 
% scaled properly. 

For q_index = l:q_max+l 

temp(:,q_index) = temp(:,q_index) / divisors(q_index); 
end 

% If combination of current v and any q provides a greater value 
% than the previous max, then remember m, Q, & V. 
if max(max(temp))>x 
X = max(max(temp)); 

[m Q] = find(temp == max(max(temp))); 

% Must do this since q starts at 0, but Matlab doesn't allow for 
% zero indexing. 

Q = Q - 1; 

V = v_values(v); 

end 

end 

% Coarse estimate of TDOA (in # of samples) 

TDOA_Coarse = Q * N1 + (-N2+1 + m) ; 

% Coarse estimate of FDOA (in Freq Bin #) 

FDOA_Coarse = V/N1*N; 

% The following 3 lines can be used to display the coarse estimates, 

% if desired. 

%disp(['The coarse TDOA estimate is: ', num2str (TDOA_Coarse) ,... 

% ' samples .']) ; 

%disp(['The coarse FDOA estimate is: ', num2str (FDOA_Coarse/N) ,... 

% ' (digital frequency).']); 

9 - kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
0 

% FINE MODE computations. 

9 - kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 
0 

S2 = conj(S2); % S2 is conjugated in basic CAE definition 
% Vector of freq "bins" to use (DON'T have to be integers!!) 
k_val = FDOA_Coarse-10 : FDOA_Coarse+10; 
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% Vectors of TDOAs to use (must be integers) 
tau_val = TDOA_Coarse-10 : TDOA_Coarse+10 ; 
done = 0; 
multiple = 1; 
decimal = 0; 

while ~done % Fine mode iterations continue until user is done. 

% Initialize to make later computations faster 
amb(length(k_val),length(tau_val))=0; 

Ntemp = N * multiple; 

for k = 1:length(k_val) % Must loop through all values of k 
% Vector of complex exponentials that will be used 
exponents = exp(-j*2*pi*k_val(k)/Ntemp*(0 :Ntemp-l)') ; 

% Must loop through all potential TDOAs 
for t = 1:length(tau_val) 

% S2 is shifted "tau" samples 
S2temp = shiftud(S2,tau_val(t),0); 

% Definition of CAF summation 

temp = abs(sum(SI.*S2temp.*exponents)) ; 

% Save CAF magnitude for the values of k & t 

amb(k,t)=temp; 

end 

end 

[k, t]=find(amb==max(max(amb))); % Find the peak of the CAF matrix 
TDOA = tau_val(t); % TDOA and FDOA associated with the peak of the 
FDOA = k_val(k); % CAF plane. These represent the final TDOA 
% & FDOA estimates. 

% The results are displayed. 

Disp ( ' ');disp(' ');disp(' '); 

disp(['The TDOA is num2str(TDOA/multiple), ' samples']); 

disp ( [ ' or num2str(TDOA/(multiple*!s)) , ' seconds.']); 

disp ( ' ') ; 

disp (['The resolution is num2str (0.5/... 

(multiple*!s)) , ' seconds.']) ; 
disp ( ' ');disp ( ' '); 

disp (['The FDOA is ', num2str (FDOA/N) , 

' in digital frequency (k/N)']); 

disp ( [ ' or ', num2str(FDOA/N*fs), ' Hz.']); disp( ' '); 

disp (['The resolution is ', num2str (0.5*... 

(lO^decimal)/N*fs), ' Hz.']); 

disp ( ' ');disp(' ');disp(' '); 

% If the signal length exceeds 524288 elements, max processing 
% capability has been achieved, and the user will not be given 
% the option of refining TDOA any further. 

If Ntemp >= 2^19 

disp('Maximum TDOA processing capability has been achieved.') 
done! = 1; 
else done! = 0; 
end 

% User chooses whether to compute more accurate TDOA &/or 
% FDOA, or to stop fine mode computations. 

Disp('Do you desire a solution with finer resolution?'); 
disp ('Select one of the following:'); disp(' '); 

if ~doneT 

disp('l. Finer resolution for TDOA.'); 

else disp(' '); 

end 
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disp('2. Finer resolution for FDOA.'); 
if ~doneT 

disp('3. Finer resolution for both TDOA and FDOA.'); 

else disp(' '); 

end 

disp ('4. The TDOA and FDOA resolutions are fine enough.'); 
disp ( ' ') ; 

choice = input('What is your selection? '); 
switch choice 

% TDOA is refined by resampling the signals at twice the 
% previous sampling rate. Increases resolution two-fold. 

Case 1 
if ~doneT 

multiple = multiple*2; 

51 = interp(Sl, 2); 

52 = interp(S2, 2); 

tau_val = TD0A*2 - 1 : TD0A*2 + 1; 

else done = 1; 

end 

clc; 

% FDOA resolution is improved by a factor of 10. 

Case 2 

decimal = decimal - 1; 

k_val = FDOA - 5*10^decimal : lO^decimal : FDOA + 5*10^decimal; 
clc; 

% Both TDOA and FDOA resolutions are improved. 

Case 3 
if ~doneT 

multiple = multiple*2; 

51 = interp(Sl, 2); 

52 = interp(S2, 2); 

tau_val = TDOA*2 - 1 : TDOA*2 + 1; 
decimal = decimal - 1; 

k_val = FDOA - 5*10^decimal : lO'^decimal : FDOA + ... 

5*10^decimal; 
else done = 1; 
end 
clc; 

otherwise 
done = 1; 
end 

if done 

disp ( ' ');disp(' '); disp('TDOA & FDOA estimation complete.'); 

end 

end 

% If user wants to see the OAF surface graphically, a call to 
% CAF_peak is made. 

Disp ( ' ');%disp(' ');disp(' '); 

choice = input... 

('Would you like to see the OAF peak graphically (Y or N)? ','s'); 
choice = upper(choice); 
switch choice 
case 'Y' 

caf_peak (Sl_orig, S2_orig, floor (TDOA/multiple) - 50, ... 
floor(TDOA/multiple) + 50, (FDOA-20)/N, (FDOA+20)/N, fs); 
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end 

TDOA = TDOA/(multiple*fs); % Returns TDOA in seconds. 
FDOA = FDOA/N*fs; % Returns FDOA in Hertz. 

Disp( 'Program Complete.'); 
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APPENDIX D. CAF II 


function [TDOA, FDOA, MaxAmb, Amb] = ... 

CAF_peak(Sl, S2, Tau_Lo, Tau_Hi, Freq_Lo, Freq_Hi, fs); 

2 - -k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

% CAF_peak(Sl, S2, Tau_Lo, Tau_Hi, Freq_Lo, Freq_Hi) takes as input: 

% two signals (SI, S2) that are row or column vectors; a range of 
% time delays (in samples) to search (Tau_Lo, Tau_Hi must be 
% integers between -N & +N); a range of digital frequencies (in 
% fractions of sampling frequency) to search (Freq_Lo, Freq_Hi must 
% be between -1/2 and 1/2, or -(N/2)/N and (N/2)/N, where N is the 
% length of the longer of the two signal vectors); and the sampling 
% frequency, fs. 

"6 

% The function computes the Cross Ambiguity Function of the two 
% signals. Four plots are produced which represent four different 
% views of the Cross Ambiguity Function magnitude versus the input 
% Tau and Frequency Offset ranges. 

"6 

% The function returns the scalars TDOA, FDOA, and MaxAmb, where 
% TDOA & FDOA are the values of Time Delay and Frequency Offset 
% that cause the Cross Ambiguity Function to peak at a magnitude 
% of MaxAmb. Amb is the matrix of values representing the CAF 
% surface. 

% Written by: LCDR Joe J. Johnson, USN 
% Last modified: 26 August 2001 

% Modified by CDR Mickey Batson, USN 
% Last modified: 10 October 2006 

9 - -k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 
0 

% Ensures that the user enters all SIX required arguments, 
if (nargin < 6) 
error... 

('6 arguments required: SI, S2, Tau_Lo, Tau_Hi, Freq_Lo, Freq_Hi'); 
end 

% Ensures that both SI & S2 are row- or column-wise vectors, 
if ((size(Sl,l)~=l)&(size(Sl,2)~=l)) | ((size(S2,1)~=1)&... 

(size(S2,2)~=1)) 

error ('SI and S2 must be row or column vectors.'); 
end 

N1 = length (SI); 

N2 = length(S2); 

51 = reshape(SI,N1,1); % SI & S2 are reshaped into column-wise 

52 = reshape(S2,N2,1); % vectors since MATLAB is more efficient 
% when manipulating columns. 

51 = [SI;zeros(N2-N1,1)]; % Ensure that SI & S2 are the same size, 

52 = [S2;zeros(N1-N2,1)]; % padding the smaller one w/ Os as neeeded. 

% This WHILE loop simply ensures that the length of SI & S2 is a power 
% of two. If not, the vectors are padded with Os until their length 
% is a power of two. This is not required, but it takes advantage of 
% the fact that MATLAB's FFT computation is significantly faster for 
% lengths which are powers of two! 

while log(length(SI))/log(2) ~= round(log(length(SI))/log(2)) 

SI(length(SI)+1) = 0; 
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S2(length(S2)+1) = 0; 
end 

N = length(SI); 

% Ensures that the Tau values entered are in the valid range, 
if abs(Tau_Lo)>N | abs(Tau_Hi)>N 

error('Tau_Lo and Tau_Hi must be in the range -N to +N. ') ; 
end 

% Ensures that Tau values entered by the user are integers, 
if (Tau_Lo ~= round(Tau_Lo)) | (Tau_Hi ~= round(Tau_Hi)) 

error('Tau_Lo and Tau_Hi must be integers.') 
end 

% Ensures that the Frequency values entered are in the valid range, 
if abs(Freq_Lo)>1/2 | abs(Freq_Hi)>1/2 

error('Freq_Lo and Freq_Hi must be in the range -.5 to +.5'); 
end 

% Ensures that the lower bounds are less than the upper bounds. 

if (Tau_Lo > Tau_Hi) | (Freq_Lo > Freq_Hi) 

error('Lower bounds must be less than upper bounds.') 

end 

% Freq values converted into integers for processing. 

Freq_Lo = round(Freq_Lo*N); 

Freq_Hi = round(Freq_Hi*N); 

% Creates vectors for the Tau & Freq values entered by the user. Used 
% for plotting... 

TauValues = [Tau_Lo:Tau_Hi]; 

FreqValues = [Freq_Lo:Freq_Hi]/N; 

% The IF statement calculates the indices required to isolate the 
% user-defined frequencies from the EFT calculations below, 
if Freq_Lo < 0 & Freq_Hi < 0 
Neg_Freq = (N+Freq_Lo+l:N+Freq_Hi+l); 

Pos_Freq = []; 

elseif Freq_Lo < 0 & Freq_Hi >= 0 
Neg_Freq = (N+Freq_Lo+l:N); 

Pos_Freq = (1:Freq_Hi+l); 
else 

Neg_Freq = [ ] ; 

Pos_Freq = (Freq_Lo+l:Freq_Hi+l); 
end 

% This FOR loop actually calculates the Cross Ambiguity Function for 
% the given range of Taus and Frequencies. Note that an FFT is 
% performed for each Tau value and then the frequencies of interest 
% are isolated using the Neg_Freq and Pos_Freq vectors obtained above. 

% For each value of Tau, the vector S2 is shifted Tau samples using a 
% call to the separate function "SHIFTUD". Samples shifted out are 
% deleted and zeros fill in on the opposite end. 

% Initializing Amb with Os makes computations much faster. 

Amb=zeros(length(Neg_Freq)tlength(Pos_Freq),length(TauValues)) ; 
for t = 1:length(TauValues) 

temp = fft ( (Si) .*conj(shiftud(S2,TauValues(t),0))); 

Amb(:,t) = [temp(Neg_Freq);temp(Pos_Freq)]; 
end 

% Only interested in the Magnitude of the Cross Ambiguity Function. 

Amb = abs(Amb); 

% The following will remove any spike that occurs at Tau = FreqOff = 0. 
% This may be desired in some cases, especially when the spike at (0,0) 
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% is due to correlation of the two signals' noise components. The 
% spike, of course, could also indicate that the two signals have no 
% TDOA or FDOA between them. 

% if find(TauValues == 0) & find(FreqValues == 0) 

% Amb(find(FreqValues==0),find(TauValues==0)) = 0; 

% end 

%clc; %Clears the MATLAB command window. 

% The four different views of the Cross Ambiguity Function plots are 
% created here. 

figure % This one is the 3-D view 

mesh(TauValues/fs,FreqValues*fs,Amb); 

xlabel('TDOA (Seconds)');ylabel('FDOA (Hertz)'); 

zlabel('Magnitude'); 

title('Cross Ambiguity Function'); 

figure 

subplot (2,1,1) % This one is the 2-D view along the TDOA axis 

mesh(TauValues/fs,FreqValues*fs, Amb) ; 

xlabel('TDOA (Seconds) ') ; 

zlabel('Magnitude' ) ; 

view(0,0) ; 

subplot (2,1,2) % This one is the 2-D view along the FDOA axis 

mesh(TauValues/fs,FreqValues*fs,Amb); 

ylabel('FDOA (Hertz) ') ; 

zlabel('Magnitude' ) ; 

view(90,0) ; 

%This one is a 2-D view looking down on the plane 
figure 

mesh(TauValues/fs,FreqValues*fs,Amb); 

xlabel('TDOA (Seconds)');ylabel('FDOA (Hertz)'); 

zlabel('Magnitude'); 

title('Cross Ambiguity Function'); 

view(0,90); 

% Finds the indices of the peak value. 

[DFO, DTO] = find(Amb==max(max(Amb))); 

TDOA = TauValues(DTO); % Finds the actual value of the TDOA. 

FDOA = FreqValues(DFO); % Finds the actual value of the FDOA. 

MaxAmb = max(max(Amb)); % Finds the actual Magnitude of the peak. 

% The remaining lines will display the numerical results of the 
% TDOA & FDOA, if desired. Since the FFT method was used for the 
% calculations, the TDOA is accurate only to within +/- 0.5 samples, 

% and the FDOA is accurate to within +/- 0.5/N in digital frequency. 

% disp ( ' '); disp ( ' ' ); 

% disp(['The TIME LAG (TDOA) is: ' , num2str(TDOA), ' Samples.']); 

% disp ( ' ' ) ; 

% disp(['The FREQ OFFSET (FDOA) is: ',num2str(FDOA),... 

% ' (Fraction of Fs).']); 

% disp(' '); disp(['Maximum Magnitude = ',num2str(MaxAmb)]); 

% disp ( ' ');disp(' - '); 

% disp('NOTE: If the CAF plot has secondary peaks whose magnitudes'); 
% disp(' are within about 80% of the Main Peak''s magnitude,'); 

% disp(' then the above results may be unreliable. Likely'); 

% disp(' reasons: The true peak is not within the range of,'); 

% disp(' Taus & Freq Offsets that you entered or the signals'); 

% disp(' may be too noisy to detect the peak.'); 
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APPENDIX E. TDOA MAIN 


% Filename: Random_MAIN.m 



% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and +posterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 

clear 
global c 

Na=0; %input('Enter number of random array samples to be generated 

: '); 

M=2; %input('Enter number of elements along x-axis, M (must be >1) 

: '); 

N=M; 
c=3e8; 
f=157e6; 
sigm = le-9; 
k = 5; 

kk = M*N - 1; 
z = [5e3;5e3]; 


wavelength=c/f; %input('Enter design wavelength of array in meters 

: '); 

L=50; %input('Enter length of sensor cluster square area in 

meters :'); 

error= 0.0; %input('Enter position error in meters :'); 

fail=0.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

p=l; 

%load('Array'); 

[x_out,x_err,y_out,y_err,distance, midpointX, 
midpointY]=Loomis_LOC_DOA(M,N,L,error); 

%distance_lamda = distance./wavelength; 
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figure (1), plot(x_err, y_err, '*r'), grid on, hold on 

vz = [0.0;0.0]; 
dt = le-6; 
zi = [2e3;le3]; 
for i = 1:9 

zco(:,i) = [midpointX (i); midpointY (i)]; 
end 

%D = reshape(D,M*N,1); 

%[peaks, iocs] = pkpick(D,Th,(M*N-accuracy)); 

%ZZC1 = zeros(2,k);%(M*N-accuracy)) ; 

%ZZC2 = zeros(2,k);%(M*N-accuracy)) ; 

%MMM = zeros(1,(M*N-accuracy)); 

%for i = 1:length(peaks) 

%peaks = sort(peaks); 

[m, zzcl, zzc2] = tdoagen(sigm, k, z, zco(:,2), distance(:,2), vz, dt); 
%ZZC1 ( :,p) = zzcl; 

%ZZC2(:,p) = zzc2; 

%MMM(:,p) = m; 

%p=p+l; 

%end 

tic 

[mi, zz, P] = nrtdoa(k,m',zzcl,zzc2,vz,sigm,zi); 

time_elapsed = toe 

%azimuth = atan2(zz (2),zz (1)); 

%Azimuth = azimuth*(180/pi) 

[m2, zzcl2, zzc22] = tdoagen(sigm, k, z, zco(:,3), distance(:,3), vz, 
dt) ; 

%compute and plot TDOA isochron contours 
%lifted from geotdof2 

%First set up extent of plot to include all features of interest 
[mi2, zz2, P2] = nrtdoa(k,m2',zzcl2,zzc22,vz,sigm,zi); 

[m3, zzcl3, zzc23] = tdoagen(sigm, k, z, zco(:,4), distance(:, 4) , vz, 

dt) ; 

[mi3, zz3, P3] = nrtdoa(k,m3',zzcl3,zzc23,vz,sigm,zi); 

[m4, zzcl4, zzc24] = tdoagen(sigm, k, z, zco(:,5), distance(:,5), vz, 

dt) ; 

[mi4, zz4, P4] = nrtdoa(k,m4',zzcl4,zzc24,vz,sigm,zi); 

[m5, zzclS, zzc25] = tdoagen(sigm, k, z, zco(:,6), distance(:, 6) , vz, 

dt) ; 

[mis, zz5, PS] = nrtdoa(k,mS',zzclS,zzc2S,vz,sigm,zi); 

[m6, zzcl6, zzc26] = tdoagen(sigm, k, z, zco(:,7), distance(:,7), vz, 

dt) ; 

[mi6, zz6, P6] = nrtdoa(k,m6',zzcl6,zzc26,vz,sigm,zi); 

[m7, zzcl7, zzc27] = tdoagen(sigm, k, z, zco(:,8), distance(:,8), vz, 

dt) ; 

[mi7, zz7, P7] = nrtdoa(k,m7',zzcl7,zzc27,vz,sigm,zi); 

[x,y]=meshgrid(-8e3:20:8e3,-8e3:20:8e3) ; 

%[Xf y]=meshgrid(-20e3:200:20e3,-20e3:200:20e3); 
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%figure (2), hold on, grid on 


xlabel('x position in meters'), ylabel('y position in meters') 
title('Line of Bearing') 

%gtext('Sensor Grid 2500 square meters: error = 0.0 meter { 

failure = 0 nodes') 

for j=l:k 

% jth contour 

tb= (1/c) . * ( ( (x-(zzcl (1, j) ) ) . "^2 + (y-(zzcl (2, j) ) ) . ^2) . ^ . 5) - ... 

((x-(zzc2 (1,j) ) ) . ^2 + (y-(zzc2 (2,j) ) ) . ^2) . ^(.5)); 
MT=max(max(tb)); 
mt=min(min(tb)); 
ttb=mi(j,1); 
tlb=[1 ttb]; 

Tb=contour(-8e3:20:8e3,-8e3:20:8e3,tb,tlb, 'b-'); 

%Tb=contour(-20e3:200:20e3,-20e3:200:20e3,tb,tlb,'b-'); 
ttb=mi(j,1)+sigm; 
tlb=[1 ttb]; 

%Tb=contour(-8e3:20:8e3,-8e3:20:8e3,tb,tlb, 'c-') ; 

%Tb=contour(-20e3:200:20e3,-20e3:200:20e3,tb,tlb,'c-'); 
ttb=mi(j,1)-sigm; 
tlb=[1 ttb]; 

%Tb=contour(-8e3:20:8e3,-8e3:20:8e3,tb,tlb, 'c-') ; 

%Tb=contour(-20e3:200:20e3,-20e3:200:20e3,tb,tlb,'c-'); 

tb2=(l/c).*(((x-(zzcl2(1,j))).^2 + (y-(zzcl2(2,j))).^2)(.5) - 
((x-(zzc22 (1,j) ) ) . ^2 + (y-(zzc22 (2,j) ) ) . ^2) . ^(.5)); 
MT2=max(max(tb2)); 
mt2=min(min(tb2)); 
ttb2=mi2(j,1); 
tlb2=[l ttb2]; 

Tb2=contour(-8e3:20:8e3,-8e3:20:8e3,tb2,tlb2, 'b-') ; 

%Tb2=contour(-20e3:200:20e3,-20e3:200:20e3,tb2,tlb2,'b-'); 
ttb2=mi2(j,1)+sigm; 
tlb2=[l ttb2]; 

%Tb2=contour(-8e3:20:8e3,-8e3:20:8e3,tb2,tlb2,'c-'); 
%Tb2=contour(-20e3:200:20e3,-20e3:200:20e3,tb,tlb, 'c -'); 
tb2=mi2(j,1)-sigm; 
tlb2=[l ttb2]; 

%Tb2=contour(-8e3:20:8e3,-8e3:20:8e3,tb2,tlb2,'c-'); 
%Tb2=contour(-20e3:200:20e3,-20e3:200:20e3,tb,tlb, 'c -'); 
tb3=(l/c).*(((X-(zzcl3(1,j))).^2 + (y-(zzcl3(2,j))).^2).^(.5) - 
((x-(zzc23 (1,j))) . ^2 + (y-(zzc23 (2,j))) . ^2) . ^ ( . 5)); 
MT3=max(max(tb3)); 
mt3=min(min(tb3)); 
ttb3=mi3(j, 1) ; 
tlb3=[1 ttb3]; 

Tb3=contour(-8e3:20:8e3,-8e3:20:8e3,tb3,tlb3, 'b-') ; 

%Tb3=contour(-20e3:200:20e3,-20e3:200:20e3,tb3,tlb3,'b-'); 
ttb3=mi3(j,1)+sigm; 
tlb3=[1 ttb3] ; 

%Tb3=contour(-8e3:20:8e3,-8e3:20:8e3,tb3,tlb3,'c-'); 
%Tb3=contour(-20e3:200:20e3,-20e3:200:20e3,tb3,tlb3, 'c-') ; 
ttb3=mi3(j,1)-sigm; 
tlb3=[1 ttb3]; 


node 


219 



%Tb3=contour(-8e3:20:8e3,-8e3:20:8e3,tb3,tlb3, 'c-'); 
%Tb3=contour(-20e3:200:20e3,-20e3:200:20e3,tb3,tlb3,'c-'); 

tb4=(l/c) .*(((x-(zzcl4(1,j))) .-2 + (y-(zzcl4(2,j))) .-2) (.5) 

((x-(zzc24 (1,j)) ) . ^2 + (y-(zzc24 (2,j) ) ) . ^2) . ^ .5)); 
MT4=max(max(tb4)); 
mt4=min(min(tb4)); 
ttb4=mi4 ( j,1); 
tlb4=[l ttb4]; 

Tb4=contour(-8e3:20:8e3,-8e3:20:8e3,tb4,tlb4,'b-'); 
%Tb4=contour(-20e3:200:20e3,-20e3:200:20e3,tb4,tlb4, 'b-'); 
ttb4=mi4(j,1)+sigm; 
tlb4=[l ttb4]; 

%Tb4=contour(-8e3:20:8e3,-8e3:20:8e3,tb4,tlb4, 'c-') ; 
%Tb4=contour(-20e3:200:20e3,-20e3:200:20e3,tb4,tlb4, 'c-'); 
ttb4=mi4(j,1)-sigm; 
tlb4=[l ttb4]; 

%Tb4=contour(-8e3:20:8e3,-8e3:20:8e3,tb4,tlb4,'c-'); 
%Tb4=contour(-20e3:200:20e3,-20e3:200:20e3,tb4,tlb4, 'c-'); 

tb5=(l/c) .* ( ( (X-(zzclS (1, j) ) ) ."-2 + (y-(zzclS (2, j) ) ) .^2) .^ ( .5) 

((x-(zzc25 (1,j))) . ^2 + (y-(zzc25 (2,j))) . ^2) . ^ ( . 5)); 
MT5=max(max(tbS)); 
mt5=min(min(tbS)); 
ttb5=mi5 ( j,1); 
tlb5=[1 ttbS]; 

Tb5=contour(-8e3:20:8e3,-8e3:20:8e3,tb5,tlb5, 'b-'); 
%Tb4=contour(-20e3:200:20e3,-20e3:200:20e3,tb4,tlb4,'b-'); 
ttb5=mi5(j,1)+sigm; 
tlb5=[1 ttbS]; 

%Tb5=contour(-8e3:20:8e3,-8e3:20:8e3,tb5,tlb5,'c-'); 
%Tb5=contour(-20e3:200:20e3,-20e3:200:20e3,tb5,tlb5, 'c-'); 
ttb5=mi5(j,1)-sigm; 
tlb5=[1 ttbS]; 

%Tb5=contour(-8e3:20:8e3,-8e3:20:8e3,tb5,tlb5,'c-'); 
%Tb5=contour(-20e3:200:20e3,-20e3:200:20e3,tb5,tlb5, 'c-'); 

tb6=(1/c) .*(((x-(zzcl6 (1,j))) .^2 + (y-(zzcl6 (2,j))) . ^2) . ^ ( . 5) 

((x-(zzc26 (1,j))) . ^2 + (y-(zzc26 (2,j))) . ^2) . ^ ( . 5)); 
MT6=max(max(tb6)); 
mt6=min(min(tb6)); 
ttb6=mi6(j,1); 
tlb6=[1 ttb6]; 

Tb6=contour(-8e3:20:8e3,-8e3:20:8e3,tb6,tlb6,'b-'); 
%Tb6=contour(-20e3:200:20e3,-20e3:200:20e3,tb6,tlb6, 'b-'); 
ttb6=mi6(j,1)+sigm; 
tlb6=[1 ttb6]; 

%Tb6=contour(-8e3:20:8e3,-8e3:20:8e3,tb6,tlb6, 'c-') ; 
%Tb6=contour(-20e3:200:20e3,-20e3:200:20e3,tb6,tlb6,'c-'); 
ttb6=mi6(j,1)-sigm; 
tlb6=[1 ttb6]; 

%Tb6=contour(-8e3:20:8e3,-8e3:20:8e3,tb6,tlb6, 'c-'); 
%Tb6=contour(-20e3:200:20e3,-20e3:200:20e3,tb6,tlb6,'c-'); 

tb7=(l/c) .*(((x-(zzcl7 (1,j))) ."2 + (y-(zzcl7(2,j))) ."2) (.5) 

((x-(zzc27 (1,j))) ."2 + (y-(zzc27 (2,j))) ."2) (.5)); 

MT7=max(max(tb7)); 
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mt7=min(min(tb7)); 
ttb7=mi7(j,1); 
tlb7=[l ttb7]; 

Tb7=contour(-8e3:20:8e3,-8e3:20:8e3,tb7,tlb7,'b-'); 
%Tb7=contour(-20e3:200:20e3,-20e3:200:20e3,tb7,tlb7, 'b-') ; 
ttb7=mi7(j,1)+sigm; 
tlb7=[l ttb7]; 

%Tb7=contour(-8e3:20:8e3,-8e3:20:8e3,tb7,tlb7,'c-'); 
%Tb7=contour(-20e3:200:20e3,-20e3:200:20e3,tb7,tlb7, 'c-') ; 
ttb7=mi7(j,1)-sigm; 
tlb7=[l ttb7]; 

%Tb7=contour(-8e3:20:8e3,-8e3:20:8e3,tb7,tlb7,'c-'); 
%Tb7=contour(-20e3:200:20e3,-20e3:200:20e3,tb7,tlb7, 'c-') ; 

End 
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APPENDIX F. TDOAGEN 


function [tdoa,zzcl,zzc2] = tdoagen(sigm,k,z,zcO,D,vz,dt) 

% Generates k component vector of tdoas 
% sigm - Standard deviation of time measurement 
% k - number of measurements 
% z - True position of emitter 
% D - Distance separating collectors 
% vz - Vector collector velocity 
% dt - Time between TDOA measurements 

% tdoa - k-component vector of simulated TDOA measurements 
% Written by HH Loomis on 7/7/99 
% Based on nrgeogrf.m version 3.41 

% s - half separation 

% zz - Matrix of collector positions, 2xk 
global c 

disp ( ' Generating') 
disp(k) 

disp ( ' TDOA Observations.') 
s=D/2; 

zz=zeros(2,k); 
for j=l:k 

%rl(:,j)=z-(zz(:,j)-[s;0]); 

%r2(:,j)=z-(zz(:,j) + [s;0]); 

%TDOA(j,1)=(1./c)*(norm(rl(:,j))-norm(r2(:,j))) + sigm*randn; 

zz(:,j)=zcO; %+vz*(j-1)*dt; 
zzcl(:,j)=zz(:,j)-zz(:,j); %[s;0]; 
zzc2(:,j)=zz(:,j)+zz(:,j) %[s;0]; 

tdoa(j,l)=(l./c)*(norm(z-(zzcl(:,j))) 

(zzc2 ( :,j)))) + sigm*randn; 

end 


norm(z- 
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APPENDIX G. NEWTON-RAPHSON 


function [mi, zz, P] = nrtdoa(k,m,zzcl,zzc2,vz,sigm,zi) 

%function computes position estimate, zi, and covariance matrix, P 
% vz - Vector collector velocity 
% zzcl - Matrix of collector 1 positions, 2xk 
% zzc2 - Matrix of collector 2 positions, 2xk 
% sigm - Standard deviation of time measurement 
% zi - iterated estimated position of emitter 
% Written by HH Loomis on 7/7/99 
% Based on nrgeogrf.m version 3.41 

% W - Conditioning matrix 
% mi - vector of TDOAs for current zi 
% A - kx2 PARTIAL DERIV. MATRIX, 

global c 

% calculate Weighting Matrix, W 

% limit=l.4*max(z); 

% [x,y]=meshgrid(-limit:.02*limit:limit,-limit:.02*limit:limit) 

W=zeros(k,k); 
for j j = l:k 

W( j j, j j) =1/ (sigm"'2) ; 

end 

% disp('W=' ) 

% disp(W) 

% start estimation loop with ii = 0 

ii=0; 

while ii<10 

% disp('Begin Iteration #:') 

% disp(ii) 


%calculate mi as function of zi, initial trial value of position 
for j j = l:k 

mi(jj,l) = (l./c)*(norm(zi-(zzcl(:,jj))) - norm(zi-(zzc2(:, jj)))) ; 
end 

% Compute rms error in mi 

dm=m'-mi; 

rms_m = norm(dm); 

% disp('mi, dm, rms_m') 

% disp(mi) 



% disp(dm) 

% disp(rms_m) 

% disp('mmi') 

% disp(mmi) 

if rms_m <= .l*sigm break 
else 

% CALCULATE PARTIAL DERIV. MATRIX, A(K,2) 

for j j = l:k 

A(jj, :) = (l/c)*(((zi-(zzcl(:,jj))) '/norm(zi ... 
-(zzcl (:,jj)))) - ... 

((zi-(zzc2(:,jj))) '/norm(zi-(zzc2 ( :,j j))))) ; 

end 

% disp('Matrix A') 

% disp(A) 

%disp ('Condition of At*W*A') 

%cond(A'*W*A) 

% COMPUTE increment in z, dzi 

P=inv(A'*W*A); 
dzi=(P)*(A'*W)*dm; 

% disp('dzi') 

% disp(dzi) 

% COMPUTE next estimate for zi 

zipl=zi+dzi; 

% disp('zipl') 

% disp(zipl) 

%ii used as loop index since i is sqrt(-l) in matlab 

% disp('end of iteration loop #') 

% disp(ii) 

ii=ii+l; 
zi=zipl; 

end %for if — else 

end %while ii 
zz=zi 

%end iteration loop 

disp('ii') 
disp(ii) 
disp('zi, dzi') 
disp(zi) 
disp(dzi) 
disp('dm') 
disp(dm) 
disp('rms_m') 
disp(rms_m) 
disp ( ' P ' ) 
disp (P) 
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APPENDIX H. RANDOM LOG DOA 


% Filename: Random_LOC.m 



% Input parameters: 

% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% L = length of square planar area (in meters) 

% error = position errors (assumed uniformly distributed between [- 
posterror and +posterror]) 


% Output: 

% X = x-coordinates of sensor nodes (in meters) 

% x_err = x-coordinates of sensor nodes with errors (in meters) 

% y = y-coordinates of sensor nodes (in meters) 

% y_err = y-coordinates of sensor nodes with errors (in meters) 

% distance = distance between the sensor nodes and the reference node 
in 

% meters 


function[x_loc,x_err,y_loc,y_err,distance,midpointX,midpointY] = 

Random_LOC(M,N,L,error,lamda) 

%M=3; 

%N=M; 

%L=10; 

%error=0; 

%X = lamda.* [0 1 2 3 4; 5 6 7 8 9; 10 11 12 13 14; 15 16 17 18 19; 20 
21 22 23 24] ; 

%Y = lamda.*[0 1234; 01234; 01234; 01234; 01234]; 

%Y = zeros(M,N); 

x_loc=[0; -50; 30; 0; -30; 25; -40; -100; 100]; 
y_loc=[0; 0; -50; 70; -50; -25; 70; 50; -80]; 

%x_loc=L*X; 

%y_loc=L*Y; 

%x_loc(1, l)=zeros(1, 1) ; 

%y_loc(1, l)=zeros(1, 1) ; 

%dx=-error+(2*error)*rand(M,N); 

%dy=-error+(2*error)*rand(M, N); 
x_err=x_loc; %+dx; 
y_err=y_loc; %+dy; 

for i = 1:9 

distance (i) = sqrt((x_err(1)-x_err(i))^2 + (y_err(1)-y_err(i))^2); 
midpointX(i) = (x_err(1)+x_err(i))/2; 
midpointY(i) = (y_err(1)+y_err(i))/2; 
end 

x=x_loc; 
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y=y_loc; 

save ARRAY x y x_err y_err distance midpointX midpointY 


% Filename: Random_LOC_DOA.m 



% Input parameters: 

% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% L = length of square planar area (in meters) 

% error = position errors (assumed uniformly distributed between [- 
posterror and +posterror]) 


% Output: 

% X = x-coordinates of sensor nodes (in meters) 

% x_err = x-coordinates of sensor nodes with errors (in meters) 

% y = y-coordinates of sensor nodes (in meters) 

% y_err = y-coordinates of sensor nodes with errors (in meters) 

% distance = distance between the sensor nodes and the reference node 
in 

% meters 


function[x_loc,x_err,y_loc,y_err,distance,midpointX,midpointY] 
Random_LOC_DOA(M,N,L,error,lamda) 

%M=3; 

%N=M; 

%L=10; 

%error=0; 


%x 

21 

= lamda.*[0 
22 23 24]; 

1234; 56789; 1011 

12 

13 

14; 

15 16 17 

18 19; 20 

%Y 

= lamda.*[0 

123 4; 0123 4; 0123 

4; 

0 1 

2 3 

4; 0 1 2 

3 4]; 


%Y = zeros(M,N); 
x_loc=L*rand(M,N)-L/2; 
y_loc=L*rand(M,N)-L/2; 

%x_loc=L*X; 

%y_loc=L*Y; 

x_loc(1, l)=zeros(1, 1) ; 
y_loc(1, l)=zeros(1, 1) ; 
dx=-error+(2*error)*rand(M,N); 
dy=-error+(2*error)*rand(M,N); 
x_err=x_loc+dx; 
y_err=y_loc+dy; 

for i = 1:M 

for j = 1:N 

distance (i, j ) = sqrt ( (x_err (1, 1)-x_err (i, j ) )'^2 + (y_err(l,l)- 

y_err(i,j))^2) ; 

midpointX (i,j) = (x_err(1,1)+x_err(i,j))/2; 
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midpointY(i,j) = (y_err(1,1)+y_err(i,j))/2; 
end 

end 

distance = distance 
midpointX = midpointX 
midpointY = midpointY 
x=x_loc; 
y=y_loc; 

save ARRAY x y x_err y_err distance midpointX midpointY 
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APPENDIX I. BEAMFORMING 


% Filename: Random_MAIN.m 



% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and +posterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 
clear 

Na=10; %input('Enter number of random array samples to be generated 

: '); 

M=5; %input('Enter number of elements along x-axis, M (must be >1) 

: '); 

N=M; 
c=3e8; 
f=157e6; 

%disp('Number of elements along y-axis, N = M') 

wavelength=c/f; %input('Enter design wavelength of array in meters 

: '); 

L=50; %input('Enter length of sensor cluster square area in 

meters :'); 

errorl= 0.0; %input('Enter position error in meters :'); 

faill=0.0; %input('Enter number of elements failed (must be less 

than MxN) : ') ; 

error2= wavelength/2; %input('Enter position error in meters :'); 
fail2=0.0; %input('Enter number of elements failed (must be less 

than MxN) : ' ) ; 

error3= wavelength/2; %input('Enter position error in meters :'); 
fail3=5.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

theta0=pi/4; %pi*rand(1,1)-pi/2; %input('Enter of AOA elevation angle 

in radians (-pi/2 < angle < pi/2) :'); 

thetal=pi/4; %theta0-pi/10 % pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) :'); 
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theta2=pi/4; %theta0+pi/10 %pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) 

phi01=pi/4 %disp('Assume AOA azimuth angle = pi/4 radians (45 degrees) 

'); 

phi02=pi/4 
phi03=pi/4 

phil=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

phi2=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

fs=500; %input('Enter sampling rate (>300) :'); 


AF_mag_sum_thetal=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta2=zeros(1,3142); 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta3=zeros(1,3142) ; 
affect the AF_cal function 

"O 

Do 

not 

change 

this. 

it 

will 


AF_mag_sum_phil=zeros(1,6284); 

AF_mag_sum_phi2=zeros(1, 6284) ; 

AF_mag_sum_phi3=zeros(1, 6284) ; 

%load('Array') ; 

%[x, x_err, y,y_err,distance,midpointX,midpointY]=Random_LOC_DOA(M,N,L,er 
ror2) ; 

% figure(l), plot(x,y, ' * ' ), grid on, hold on 
for j=l:1:Na 

j 

[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 
x(l,1)=0.0; 
x_err(1,1)=0.0; 
y(1,1)=0.0; 
y_err(1,1)=0.0; 

figure(l), plot(x,y, ' * ' ), grid on, hold on 

[mag_Wtsl,ang_Wtsl]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO1, phil, phi2, fs 
,wavelength,f); 

[mag_Wts_faill]=Random_FAIL(N,M,mag_Wtsl, faill); 

[THETAl, 

AF_mag_thetal]=AF_cal_theta(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
phiO 1) ; 

[PHIl, 

AF_mag_phil]=AF_cal_phi(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
thetaO); 

AF_mag_sum_t het a1=AF_mag_sum_t het a1+AF_mag_t het a1; 

AF_mag_sum_phil=AF_mag_sum_phi1+AF_mag_phi1; 


[mag_Wts2,ang_Wts2]=LMS(N,M,x,y,thetaO,thetal,theta2,phi02,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail2]=Random_FAIL(N,M,mag_Wts2,fail2); 

[THETA2, 

AF_mag_theta2]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,wave 
length, phi02); 

[PHI2, 

AF_mag_phi2]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,waveleng 
th, thetaO); 
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AF_mag_sum_theta2=AF_mag_sum_theta2+AF_mag_theta2; 
AF_mag_sum_phi2=AF_mag_sum_phi2+AF_mag_phi2; 


[mag_Wts3,ang_Wts3]=LMS(N,M,x,y,thetaO,thetal,theta2,phi03,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail3]=Random_FAIL(N,M,mag_Wts3,fail3); 

[THETA3, 

AF_mag_theta3]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,wave 
length, phi03); 

[PHI3, 

AF_mag_phi3]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,waveleng 
th, thetaO); 

AF_mag_sum_theta3=AF_mag_sum_theta3+AF_mag_theta3; 

AF_m a g_ s um_p h 13 = AF_m a g_ s um_p h 13 + AF_m a g_p h 13; 


end 

xlabel('x-direction in meters (m)'), ylabel('y-direction in meters 
(m) ' ) 

title('50 square meter grid') 
avg_AF_mag_thetal=AF_mag_sum_thetal; 
avg_AF_mag_phil=AF_mag_sum_phi1; 

%avg_AF_mag_theta2=AF_mag_sum_theta2; 
avg_AF_mag_phi2=AF_mag_sum_phi2; 

%avg_AF_mag_theta3=AF_mag_sum_theta3; 
avg_AF_mag_phi3=AF_mag_sum_phi3; 

theta=(THETAl/pi)*180; 
phil=(PHIl/pi)*180; 
phi2=(PHI2/pi)*180; 
phi3=(PHI3/pi)*180; 

theta_degrees=(thetaO/pi)*180 
phi_degreesl=(phiOl/pi)*180 
phi_degrees2=(phi02/pi)*180 
phi_degrees3=(phi03/pi)*180 

% 2D plot 


figure (2), subplot(3,1,1) 

plot(phil, 20*logl0(avg_AF_mag_phil),'r'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), 

xlim([0 360]),ylim([-10 30]) 
subplot (3,1,2) 

plot(phi2,20*logl0(avg_AF_mag_phi2),'b'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), 

xlim([0 360]),ylim([-10 30]) 
title('l meter error'); 
subplot(3,1,3) 

plot(phi3,20*logl0(avg_AF_mag_phi3),'g'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), 
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xlim([0 360]),ylim([-10 30]) 

title]'1 meter error and 20 nodes fail') 

plotedit on 

zoom on 

orient tall 

% 2D Polar plot 


figure (3), polardbl(THETAl, 20*logl0(avg_AF_mag_thetal),-60, 'm'), grid 
on 

orient tall 

figure (4), polar(PHIl,20*logl0(avg_AF_mag_phil)r'), grid on, hold on 
polar(PHI2,20*logl0(avg_AF_mag_phi2), 'b') , 
polar(PHI3, 20*logl0(avg_AF_mag_phi3), 'g') 

legend('No errors and no failures','\lambda/2 error in node location', 
'\lambda/2 error in node location and 5 node failures') 
title('Adaptive Beamforming - 157 MHz') 

gtext('50 square meter grid, where \lambda = 1.9108 m') 
orient tall 

%PG_plot(THFTA,avg_AF_mag_theta) 

%PG_plot(PHI,avg_AF_mag_phi) 

%PG_plot_polar(THETA,avg_AF_mag) 

%load('AF_256_0deg_avg.mat','THETA','avg_AF_mag'); 
polardbl(THETA,20*logl0(avg_AF_mag),-60, 'm'); 
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APPENDIX!. RANDOM LOG 


% Filename: Random_LOC.m 



% Input parameters: 

% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% L = length of square planar area (in meters) 

% error = position errors (assumed uniformly distributed between [- 
posterror and +posterror]) 


% Output: 

% X = x-coordinates of sensor nodes (in meters) 

% x_err = x-coordinates of sensor nodes with errors 
% y = y-coordinates of sensor nodes (in meters) 

% y_err = y-coordinates of sensor nodes with errors 


(in meters) 
(in meters) 


function[x,x_err,y,y_err] = Random_LOC(M,N,L,error) 

x=L*rand(M,N); 
y=L*rand(M,N); 

dx=-error+(2*error)*rand(M,N); 
dy=-error+(2*error)*rand(M,N); 
x_err=x+dx; 
y_err=y+dy; 
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APPENDIX K. 3D BEAMFORMING 


% Filename: BeamSD.m 



% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and +posterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 
clear 

Na=10; %input('Enter number of random array samples to be generated 

:' ) ; 

M=5; %input('Enter number of elements along x-axis, M (must be >1) 

:' ) ; 

N=M; 
c=3e8; 
f=157e6; 

%disp('Number of elements along y-axis, N = M') 

wavelength=c/f; %input('Enter design wavelength of array in meters 

:' ) ; 

L=10; %input('Enter length of sensor cluster square area in 

meters :'); 

H=1000; % max array height 
velocity = 10.0; 

errorl= 0.0; %input('Enter position error in meters :'); 

faill=0.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

theta0=-pi/4; %pi*rand(1,1)-pi/2; %input('Enter of AOA elevation angle 
in radians (-pi/2 < angle < pi/2) :'); 

thetal=0.0; 
theta2=0.0; 

phi01=pi/4; %disp('Assume AOA azimuth angle = pi/4 radians (45 degrees) 

'); 

phil=0.0; %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 
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phi2=0.0; 

fsl=500; %input('Enter sampling rate (>300) 

AF_mag_sum_thetal=zeros(1,3142); % Do not change this, it will 

affect the AF_cal function 
AF_mag_sum_phil=zeros(1, 6284) ; 

%load('Array3D' ) ; 

% Establishes the initial node positions 

[x,x_err,y,y_err,z,z_err]=Random_LOC_DOA3D(M,N,H,L,errorl); 

%x_err(1,1)=0.0; 

%x=x_err; 

%y_err(1,1)=0.0; 

%y=y_err; 

%z=z_err; 

tic 

for j=l:1:Na 

j 

%[x,x_err,y,y_err,z,z_err]=Random_LOC_DOA3D(M,N,H,L,errorl); 

%x (1, 1) =0.0; 

%x_err(1,1)=0.0; 

%y(1,1)=0.0; 

%y_err(1,1)=0.0; 

%z_err(1,1)=0.0; 

figure(1), plot3(x,y, z,'*'), grid on, hold on 

[mag_Wtsi,ang_Wtsl]=LMS3D(N,M,H,x,y,z,thetaO,thetal,theta2,phiOl,phil,p 
hi2,fsi,wavelength,f); 

[mag_Wts_faill]=Random_FAIL(N,M,mag_Wtsl,faill); 

[THETAl, 

AF_mag_thetal]=AF_cal_theta3D(N,M,x,y,z,mag_Wts_faill, ang_Wtsl, waveleng 
th, phiO1); 

[PHIl, 

AF_mag_phil]=AF_cal_phi3D(N,M,x,y,z,mag_Wts_faill,ang_Wtsi,wavelength, 
thetaO); 

AF_mag_sum_t het a1=AF_mag_sum_t het a1+AF_mag_t het a1; 

AF_mag_sum_phil=AF_mag_sum_phi1+AF_mag_phi1; 

[x,x_err,y,y_err,z,z_err]=Random_LOC_DOA3Dvelocity(M,N,H,x,y,z,velocity 

); 


end 

avg_AF_mag_t het a1=AF_mag_sum_t het a1; 
avg_AF_mag_phil=AF_mag_sum_phi1; 
time_elapsed = toe 

thetal=(THETAl/pi)*180; 
phil=(PHIl/pi)*180; 

% 2D Polar plot 


%figure (2), polardbl(THETAl, 20*logl0(avg_AF_mag_thetal),-60,'m'), grid 
on 

%orient tall 

figure (2), polar(THETAl,20*logl0(avg_AF_mag_thetal),'r'), grid on, 
hold on 
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axis auto 


title('Adaptive Beamforming - 3D') 

%gtext('10 square meter grid') 
orient tall 

figure (3), polar(PHIl,20*logl0(avg_AF_mag_phil),'r'), grid on, hold on 
axis auto 

title('Adaptive Beamforming - 3D') 

%gtext( '25 square meter grid') 
orient tall 
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APPENDIX L. RANDOM LOG VELOCITY 


% Filename: Random_L0C3Dvelocity.m 


% Input parameters: 

% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% L = length of square planar area (in meters) 

% error = position errors (assumed uniformly distributed between [- 
posterror and fposterror] ) 


% Output: 

% X = x-coordinates of sensor nodes (in meters) 

% x_err = x-coordinates of sensor nodes with errors (in meters) 

% y = y-coordinates of sensor nodes (in meters) 

% y_err = y-coordinates of sensor nodes with errors (in meters) 

% distance = distance between the sensor nodes and the reference node 
in 

% meters 


function[x_loc,x_err,y_loc,y_err,z_loc,z_err,distance,midpointX,midpoin 
tY, midpointZ] = Random_LOC(M,N,H,x,y,z,velocity) 

%M=3; 

%N=M; 

%L=10; 

%error=0; 


%x 

= lamda.*[0 1234; 

5 6 7 8 

9; 10 11 12 13 14; 15 16 17 18 19; 20 

21 

22 23 24] ; 



%Y 

= lamda.*[0 1 2 3 4; 0 

12 3 4 

; 01234; 01234; 01234]; 

%Y 

= zeros(M,N); 



dx= 

^-velocityf(2*velocity) 

*rand(M, 

N) ; 

dy= 

^-velocityf(2*velocity) 

*rand(M, 

N) ; 

dz = 

^-velocityf(2*velocity) 

*rand(M, 

N) ; 


x_loc=x+dx; %+rand(M,N); 
y_loc=y+dy; %+rand(M,N); 
z_loc=z+dz; %+rand(M,N); 

%x_loc=L*X; 

%y_loc=L*Y; 

x_loc(1,l)=zeros(1, 1) ; 
y_loc(1, l)=zeros(1, 1) ; 

x_err=x_loc;%+dx; 
y_err=y_loc;%+dy; 
z_err=z_loc;%+dz; 

for i = 1:M 

for j = 1:N 
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distance (i, j ) = sqrt ( {x_err (1, 1)-x_err (i, j ) ) "'2 + (y_err(l,l)- 

y_err(i,j))^2 + (z_err(1,1)-z_err(i,j))^2); 

midpointX(i,j) = (x_err(1,1)+x_err(i,j))/2; 
midpointY(i,j) = (y_err(1,1)+y_err(i,j))/2; 
midpointZ(i,j) = (z_err(1,1)+z_err(i,j))/2; 
end 

end 

distance = distance; 
midpointX = midpointX; 
midpointY = midpointY; 
midpointZ = midpointZ; 
x=x_loc; 
y=y_loc; 
z=z_loc; 

%save ARRAY3D x y z x_err y_err z_err distance midpointX midpointY 
midpointZ 
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APPENDIX M. BEAMSCAN 



% Filename: Random_MAIN.m 



% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and fposterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 
clear 

Na=10; %input('Enter number of random array samples to be generated 

:' ) ; 

M=5; %input('Enter number of elements along x-axis, M (must be >1) 

:' ) ; 

N=M; 
c=3e8; 
f=157e6; 

%disp('Number of elements along y-axis, N = M') 

wavelength=c/f; %input('Enter design wavelength of array in meters 

:' ) ; 

L=100*wavelength; %input('Enter length of sensor cluster square 

area in meters :' ) ; 


error1= 0.0; 
faill=0.0; 

%input( 'Enter position 
%input( 'Enter number 

error in meters :'); 
of elements failed 

(must 

be 

less 

than MxN) :') 
error2= 0.0; 
fail2=0.0; 

* } 

%input( 'Enter position 
%input('Enter number 

error in meters :'); 
of elements failed 

(must 

be 

less 

than MxN) :') 
error3= 0.0; 
fail3=0.0; 

* } 

%input( 'Enter position 
%input('Enter number 

error in meters :'); 
of elements failed 

(must 

be 

less 

than MxN) :') 
error4= 0.0; 
fail4=0.0; 

* } 

%input( 'Enter position 
%input('Enter number 

error in meters :'); 
of elements failed 

(must 

be 

less 

than MxN) :') 
error5= 0.0; 
fail5=0.0; 

* } 

%input( 'Enter position 
%input('Enter number 

error in meters :'); 
of elements failed 

(must 

be 

less 


than MxN) : ') ; 
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error6= 0.0; %input('Enter position error in meters 

fail6=0.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

error7= 0.0; %input('Enter position error in meters :'); 

fail7=0.0; %input('Enter number of elements failed (must be less 

than MxN) : ' ) ; 


theta0=pi/10; %pi*rand(1,1)-pi/2 ; %input('Enter of AOA elevation angle 
in radians (-pi/2 < angle < pi/2) :'); 

thetal=pi/10; %theta0-pi/10 % pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) :'); 

theta2=pi/10; %theta0+pi/10 %pi*rand(1,1)-pi/2; %input('Enter of AOA 

elevation angle in radians (-pi/2 < angle < pi/2) :'); 

theta3=pi/10; 

theta4=pi/10; 

theta5=pi/10; 

theta6=pi/10; 

phi01=pi/3; %disp('Assume AOA azimuth angle = pi/4 radians (45 degrees) 

') ; 

phi02=pi/4; 
phi03=pi/6; 
phi04=pi/7; 
phi05=pi/8; 
phiO 6=pi/9; 
phi07=pi/10; 

phil=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

phi2=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

fs=500; %input('Enter sampling rate (>300) :'); 


AF_mag_sum_thetal=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta2=zeros(1,3142); 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta3=zeros(1,3142) ; 
affect the AF_cal function 

"O 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta4=zeros(1,3142); 
affect the AF_cal function 

"O 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta5=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta6=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta7=zeros(1,3142); 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 


AF_mag_sum_phil=zeros(1,6284); 

AF_mag_sum_phi2=zeros(1, 6284) ; 

AF_mag_sum_phi3=zeros(1, 6284) ; 

AF_mag_sum_phi4=zeros(1,6284); 

AF_mag_sum_phi5=zeros(1, 6284); 

AF_mag_sum_phi6=zeros(1, 6284) ; 

AF_mag_sum_phi7=zeros(1, 6284) ; 

%load('Array' ) ; 

%[x, x_err, y, y_err, distance, midpointX, midpointY]=Random_LOC_DOA(M,N,L,er 
ror2) ; 
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% figure(1), plot(x,y,'*'), grid on, hold on 
for j=l:1:Na 

j 

[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 
x(l,1)=0.0; 
x_err(1,1)=0.0; 
y(1,1)=0.0; 
y_err(1,1)=0.0; 

figure(1), plot(x,y,'*'), grid on, hold on 
tic 

[mag_Wtsi,ang_Wtsl]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO1,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_faill]=Random_FAIL(N,M,mag_Wtsl,faill); 

[THETAl, 

AF_mag_thetal]=AF_cal_theta(N,M,x,y,mag_Wts_faill,ang_Wtsi,wavelength, 
phiO 1) ; 

[PHIl, 

AF_mag_phil]=AF_cal_phi(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
thetaO); 

AF_mag_sum_t het a1=AF_mag_sum_t het a1+AF_mag_t het a1; 

AF_mag_sum_phil=AF_mag_sum_phi1+AF_mag_phi1; 
time_elapsedl = toe 

tic 

[mag_Wts2,ang_Wts2]=LMS(N,M,x,y,thetaO,thetal,theta2,phi02,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail2]=Random_FAIL(N,M,mag_Wts2,fail2); 

[THETA2, 

AF_mag_theta2]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2, wave 
length, phi02); 

[PHI2, 

AF_mag_phi2]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,waveleng 
th, thetaO); 

AF_mag_sum_theta2=AF_mag_sum_theta2+AF_mag_theta2; 
AF_mag_sum_phi2=AF_mag_sum_phi2+AF_mag_phi2; 
time_elapsed2 = toe 

tic 

[mag_Wts3,ang_Wts3]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO3,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail3]=Random_FAIL(N,M,mag_Wts3, fail3); 

[THETA3, 

AF_mag_theta3]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,wave 
length, phi03); 

[PHI3, 

AF_mag_phi3]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,waveleng 
th, thetaO); 

AF_mag_sum_theta3=AF_mag_sum_theta3+AF_mag_theta3; 
AF_mag_sum_phi3=AF_mag_sum_phi3+AF_mag_phi3; 
time_elapsed3 = toe 

tic 
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[mag_Wts4,ang_Wts4]=LMS(N,M,x,y,thetaO,thetal,theta2,phi04,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail4]=Random_FAIL(N,M,mag_Wts4,fail4); 

[THETA4, 

AF_mag_theta4]=AF_cal_theta(N,M,x,y,mag_Wts_fail4,ang_Wts4,wavelength, 
phi04) ; 

[PHI4, 

AF_mag_phi4]=AF_cal_phi(N,M,x,y,mag_Wts_fail4,ang_Wts4,wavelength, 
thetaO); 

AF_mag_sum_theta4=AF_mag_sum_theta4+AF_mag_theta4; 
AF_mag_sum_phi4=AF_mag_sum_phi4+AF_mag_phi4; 
time_elapsed4 = toe 

tic 

[mag_Wts5,ang_Wts5]=LMS(N,M,x,y,thetaO,thetal,theta2,phiOS,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_falls]=Random_FAIL(N,M,mag_Wts5, failS); 

[THETAS, 

AF_mag_theta5]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_failS,ang_Wts5,wave 
length, phiOS); 

[PHIS, 

AF_mag_phiS]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_failS,ang_WtsS,waveleng 
th, thetaO); 

AF_mag_sum_thetaS=AF_mag_sum_thetaS+AF_mag_thetaS; 
AF_mag_sum_phiS=AF_mag_sum_phiS+AF_mag_phiS; 
time_elapsedS = toe 

tic 

[mag_Wts6,ang_Wts6]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO 6,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail6]=Random_FAIL(N,M,mag_Wts6,failG); 

[THETA6, 

AF_mag_theta6]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail6,ang_Wts6,wave 
length, phi06); 

[PHI6, 

AF_mag_phi6]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail6,ang_Wts6,waveleng 
th, thetaO); 

AF_mag_sum_theta6=AF_mag_sum_theta6+AF_mag_theta6; 
AF_mag_sum_phi6=AF_mag_sum_phi6+AF_mag_phi6; 
time_elapsed6 = toe 

tic 

[mag_Wts7,ang_Wts7]=LMS(N,M,x,y,thetaO,thetal,theta2,phi07,phil, phi2, fs 
,wavelength,f); 

[mag_Wts_fail7]=Random_FAIL(N,M,mag_Wts7, fail7); 

[THETA7, 

AF_mag_theta7]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail7,ang_Wts7,wave 
length, phi07); 

[PHI7, 

AF_mag_phi7]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail7,ang_Wts7,waveleng 
th, thetaO); 

AF_mag_sum_theta7=AF_mag_sum_theta7+AF_mag_theta7; 
AF_mag_sum_phi7=AF_mag_sum_phi7+AF_mag_phi7; 
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time_elapsed7 = toe 


end 

xlabel( 'x-direction in meters (m)'), ylabel( 'y-direction in meters 
(m) ' ) 

title ('lOOMambda square meter grid') 
avg_AF_mag_t het a 1 =AF_mag_sum_t het a 1; 
avg_AF_mag_phil=AF_mag_sum_phi1; 

%avg_AF_mag_theta2=AF_mag_sum_theta2; 
avg_AF_mag_phi2=AF_mag_sum_phi2; 

%avg_AF_mag_theta3=AF_mag_sum_theta3; 
avg_AF_mag_phi3=AF_mag_sum_phi3; 
avg_AF_mag_phi4=AF_mag_sum_phi4; 
avg_AF_mag_phi5=AF_mag_sum_phi5; 
avg_AF_mag_phi6=AF_mag_sum_phi6; 
avg_AF_mag_phi7 =AF_mag_sum_phi7; 

%fname = input( 'Save average array factor data for elevation angle as 
(with .mat extension) : 's'); 

%save(fname,'THETA','avg_AF_mag_theta') 

%fname = input( 'Save average array factor data for azimuth as (with 
.mat extension) : ', 's'); 

%save(fname,'PHI','avg_AF_mag_phi' ) 

theta=(THETAl/pi)*180; 
phil=(PHIl/pi)*180; 
phi2=(PHI2/pi)*180; 
phi3=(PHI3/pi)*180; 
phi4=(PHI4/pi)*180; 
phi5=(PHI5/pi)*180; 
phi6=(PHI6/pi)*180; 
phi7=(PHI7/pi)*180; 

theta_degrees=(thetaO/pi)*180 
phi_degreesl=(phiOl/pi)*180 
phi_degrees2=(phi02/pi)*180 
phi_degrees3=(phi03/pi)*180 
phi_degrees4=(phi04/pi)*180 
phi_degrees5=(phi05/pi)*180 
phi_degrees6=(phi06/pi)*180 
phi_degrees7=(phi07/pi)*180 

% 2D plot 


figure (2), polar(PHIl,20*logl0(avg_AF_mag_phil),'r'), grid on, hold on 

polar(PHI2,20*logl0(avg_AF_mag_phi2),'b') , 

polar(PHI3,20*logl0(avg_AF_mag_phi3) , ' g' ) , 

polar(PHI4,20*logl0(avg_AF_mag_phi4) , ' k' ) 

polar(PHIS,20*logl0(avg_AF_mag_phi5),'m'), 

polar(PHI6,20*logl0(avg_AF_mag_phi6) , ' y' ) , 

polar(PHI7,20*logl0(avg_AF_mag_phi7),'c') 

legend('\lambda/3 angle-of-arrival', '\lambda/4 angle-of-arrival', 
'\lambda/6 angle-of-arrival','\lambda/7 angle-of-arrival', '\lambda/8 
angle-of-arrival', '\lambda/9 angle-of-arrival','\lambda/10 angle-of- 
arrival' ) 

title('Adaptive Beamforming - 157 MHz') 
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gtext ('lOOMambda square meter grid, where \lambda = 1.9108 m') 
orient tall 

%PG_plot(THETA,avg_AF_mag_theta) 

%PG_plot(PHI,avg_AF_mag_phi) 

%PG_plot_polar(THETA,avg_AF_mag) 

%load('AF_256_0deg_avg.mat','THETA','avg_AF_mag'); 
polardbl(THETA,20*logl0(avg_AF_mag),-60, 'm' ); 
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APPENDIX N. LEAST MEAN SQUARE (LMS) 


"6 


% Filename: LMS.m 
"0 


% Input parameters: 

% N = number of elements along y-axis (N must be equal to M) 

% M = number of elements along x-axis (M must be equal to N) 

% X = x-coordinates of sensor nodes (in meters) 

% y = y-coordinates of sensor nodes (in meters) 

% thetaO = AOA theta angle (assume AOA phi angle = 45 degrees) 
% fs = sampling rate for the LMS algorithm 
"0 


% Output: 

% mag_Wts = magnitude of the weights 
% ang_Wts = angle of the weights (in radians) 
"0 


function[mag_Wts,ang_Wts] 

LMS(N,M,X,y,thetaO,thetal,theta2,phiO,phil,phi2,fs,wavelength,f) 

"0 

% Frequency settings 
"6 

c=3e8; 

mfreqO=c/wavelength; 

= 100 Hz 
mfreql=mfreqO; 
fl 

mfreq2=mfreqO; 
f2 

fsamp=2*mfreqO; 
lamdad=wavelength; 
lamdaO=wavelength; 
modulating signal 
lamdal=wavelength; 
modulating signal fl 
lamda2=wavelength; 
modulating signal f2 

"0 

% Steering vector for desired signal 
"6 

% theta0=0; % main beam arriving 

elevation angle for desired signal 

%phi0=pi/4; % main beam 

arriving azimuth angle for desired signal 
for nO=l:1:N 

for mO=l:1:M 


% set freq of desired modulating signal 
% set freq of interference modulating signal 
% set freq of interference modulating signal 

% wavelength of design freq of the array 


% wavelength 

of carrier 

freq 

of desired 

wavelength of 

carrier 

freq 

of 

interference 

wavelength of 

carrier 

freq 

of 

interference 
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Array_vectorO(mO,nO)=exp(j *(2*pi/lamda0)*(x(mO,nO) 
)+y(mO,nO)*sin(thetaO)*sin(phiO))); 
end 

end 

"6 

% steering vector for interference signal fl 
"6 

thetal=pi/3; 

elevation angle for interference signal fl 
%phil=pi/4; 

arriving azimuth angle for interference signal fl 
for nl=l:1:N 

for ml=l:1:M 

Array_vector1(ml,nl)=exp(j *(2*pi/lamdal)*(x(ml,nl) 
)+y(ml,nl)*sin(thetal)*sin(phil))); 
end 

end 

"6 

% steering vector for interference signal f2 
"0 

theta2=-pi/3; 

elevation angle for interference signal f2 
%phi2=pi/4; 

arriving azimuth angle for interference signal f2 
for n2=l:1:N 

for m2=l:1:M 

Array_vector2(m2,n2)=exp(j *(2*pi/lamda2)*(x(m2,n2) 
)+y(m2,n2)*sin(theta2)*sin(phi2))); 
end 

end 

"6 

% LMS algorithm 
"0 

% fs=500; 

w(N,M,l)=zeros; 

weights to zero 

conv_factor=0.001; 

for LMS algorithm 

phase0=2*pi*rand(1, fs) ; 

random phase for desired signal 

phasel=2*pi*rand(1, fs); 

random phase for interference signal fl 
phase2=2*pi*rand(1,fs); 

random phase for interference signal f2 
sigma=le-6; %0.1 * set the noise power 
for t=l:1:fs; 

noise(:,:,t)=sigma*randn(M,N); 
r (t)=cos(2*pi*mfreqO*t/fsamp+phaseO(t)) ; 
reference signal 

s(t)=cos(2*pi*mfreqO*t/fsamp+phaseO(t)) ; 
Desired modulating signal 


*sin(thetaO)*cos(phiO 


%main beam arriving 
% main beam 

*sin(thetal)*cos(phil 

% main beam arriving 
% main beam 

*sin (theta2)*cos(phi2 


% sampling rate 

% set initial 

% convergent factor 
% generate uniform 
% generate uniform 
% generate uniform 
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f1(t)=cos(2*pi*mfreql*t/fsamp+phasel (t)); 
interference modulating signal fl 

f2(t)=cos(2*pi*mfreq2*t/fsamp+phase2(t)); 
interference modulating signal f2 

X ( :, :,t)=s(t)*Array_vectorO + f1(t)*Array_vector1 

f2(t)*Array_vector2 + noiset); 

e(t)=conj(r(t))-sum(sum((conj(X(:,:,t))).*w(:,:,t))); 
calculate the conjugate of the error 

w (:, :,t + 1)=w(:, :,t) + conv_factor*X(:, :, t)*e(t) ; 

calculate new weights for next iteration 

eig_value(:,t)=eig(X(:,:,t)*(conj(X(:,:,t)).')); 
calculate eigenvalue of xx* during each sample 
max_eigvalue(t)=max(eig_value(:,t)); 
calculate maximum eigenvalue of xx* during each sample 
end 

"6 

% Graphical results display 
"6 

%square_error=abs((e(1:length(e)).*e(1:length(e)))); 

% plot(1:length(e),square_error) 

%Mean_sq_error=sum(square_error)/fs 
% Title('Convergence of the LMS algorithm') 

% xlabel('No. of samples'), ylabel('Magnitude(Error^2)'), xlim([0 fs]) 
grid on 
"6 

% Numerical results display 
"0 

%Max_convergence_factor = 2/max(max_eigvalue); 
mag_Wts=abs(conj(w(: , :,t))); 
ang_Wts=angle(conj(w(: , :,t))); 
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APPENDIX O. LEAST MEAN SQUARE 3D 


"6 


% Filename: LMS3D.m 
"0 


% Input parameters: 

% N = number of elements along y-axis (N must be equal to M) 

% M = number of elements along x-axis (M must be equal to N) 

% X = x-coordinates of sensor nodes (in meters) 

% y = y-coordinates of sensor nodes (in meters) 

% thetaO = AOA theta angle (assume AOA phi angle = 45 degrees) 
% fs = sampling rate for the LMS algorithm 
"0 


% Output: 

% mag_Wts = magnitude of the weights 
% ang_Wts = angle of the weights (in radians) 
"0 


function[mag_Wts,ang_Wts] 

LMS(N,M,H,x,y,z,thetaO,thetal,theta2,phiO,phil,phi2,fs,wavelength,f) 

"0 

% Frequency settings 
"6 

c=3e8; 

mfreqO=c/wavelength; 

= 100 Hz 
mfreql=mfreqO; 
fl 

mfreq2=mfreqO; 
f2 

fsamp=2*mfreqO; 
lamdad=wavelength; 
lamdaO=wavelength; 
modulating signal 
lamdal=wavelength; 
modulating signal fl 
lamda2=wavelength; 
modulating signal f2 

"0 

% Steering vector for desired signal 
"6 

% theta0=0; % main beam arriving 

elevation angle for desired signal 

%phi0=pi/4; % main beam 

arriving azimuth angle for desired signal 
for nO=l:1:N 

for mO=l:1:M 


% set freq of desired modulating signal 
% set freq of interference modulating signal 
% set freq of interference modulating signal 

% wavelength of design freq of the array 


% wavelength 

of carrier 

freq 

of desired 

wavelength of 

carrier 

freq 

of 

interference 

wavelength of 

carrier 

freq 

of 

interference 
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Array_vectorO(mO,nO)=exp(j *(2*pi/lamda0)*(x(mO,nO) 
)+y(mO,nO)*sin(thetaO)*sin(phiO)+z(mO,nO))); 
end 

end 

"6 

% steering vector for interference signal fl 
"6 

thetal=pi/3; 

elevation angle for interference signal fl 
%phil=pi/4; 

arriving azimuth angle for interference signal fl 
for nl=l:1:N 

for ml=l:1:M 

Array_vector1(ml,nl)=exp(j *(2*pi/lamdal)*(x(ml,nl) 
)+y(ml,nl)*sin(thetal)*sin(phil)+z(ml,nl))); 
end 

end 

"6 

% steering vector for interference signal f2 
"0 

theta2=-pi/3; 

elevation angle for interference signal f2 
%phi2=pi/4; 

arriving azimuth angle for interference signal f2 
for n2=l:1:N 

for m2=l:1:M 

Array_vector2(m2,n2)=exp(j *(2*pi/lamda2)*(x(m2,n2) 
)+y(m2,n2)*sin(theta2)*sin(phi2)+z(m2, n2))) ; 
end 

end 

"6 

% LMS algorithm 
"0 

% fs=500; 

w(N,M,l)=zeros; 

weights to zero 

conv_factor=0.001; 

for LMS algorithm 

phase0=2*pi*rand(1, fs) ; 

random phase for desired signal 

phasel=2*pi*rand(1, fs); 

random phase for interference signal fl 
phase2=2*pi*rand(1,fs); 

random phase for interference signal f2 
sigma=le-6; %0.1 * set the noise power 
for t=l:1:fs; 

noise(:,:,t)=sigma*randn(M,N); 
r (t)=cos(2*pi*mfreqO*t/fsamp+phaseO(t)) ; 
reference signal 

s(t)=cos(2*pi*mfreqO*t/fsamp+phaseO(t)) ; 
Desired modulating signal 


*sin(thetaO)*cos(phiO 


%main beam arriving 
% main beam 

*sin(thetal)*cos(phil 

% main beam arriving 
% main beam 

*sin (theta2)*cos(phi2 


% sampling rate 

% set initial 

% convergent factor 
% generate uniform 
% generate uniform 
% generate uniform 
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f1(t)=cos(2*pi*mfreql*t/fsamp+phasel(t)); 
interference modulating signal fl 

f2(t)=cos(2*pi*mfreq2*t/fsamp+phase2(t)) ; 
interference modulating signal f2 

X(:,:,t)=s(t)*Array_vectorO + f1(t)*Array_vector1 

f2(t)*Array_vector2 + noiset) ; 

e(t)=conj(r(t))-sum(sum((conj(X(:, :,t))) .*w(:, :,t))) ; 
calculate the conjugate of the error 

w(:, :,t + 1)=w(:, :,t) + conv_factor*X(:, :,t)*e (t); 

calculate new weights for next iteration 

eig_value(:,t)=eig(X(:,:,t)*(conj(X(:,:,t)).')); 
calculate eigenvalue of xx* during each sample 
max_eigvalue(t)=max(eig_value(:,t)); 
calculate maximum eigenvalue of xx* during each sample 
end 

"6 

% Graphical results display 
"6 

%square_error=abs((e(1:length(e)).*e(1:length(e)))); 

% plot(1:length(e),square_error) 

%Mean_sq_error=sum(square_error)/fs 
% Title('Convergence of the LMS algorithm') 

% xlabel('No. of samples'), ylabel('Magnitude(Error^2)') , xlim([0 fs]) 
grid on 
"6 

% Numerical results display 
"0 

%Max_convergence_factor = 2/max(max_eigvalue); 

mag_Wts=abs(conj(w(:,:,t)));ang_Wts=angle(conj(w(:,:,t))); 
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APPENDIX P. BEAMWIDTH 


% Filename: Beamwidth.m 


% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and +posterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 
clear 

Na=20; %input('Enter number of random array samples to be generated 

: '); 

M=5; %input('Enter number of elements along x-axis, M (must be >1) 

: '); 

N=M; 
c=3e8; 
f=157e6; 

%disp('Number of elements along y-axis, N = M') 

wavelength=c/f; %input('Enter design wavelength of array in meters 

: '); 

L=5*wavelength; %input('Enter length of sensor cluster square 

area in meters : ' ); 

errorl= 0.0; %input('Enter position error in meters :'); 

faill=0.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

error2= 0.0; %wavelength/2; %input('Enter position error in meters 

: '); 

fail2=0.0; %input('Enter number of elements failed (must be less 

than MxN) : ' ); 

error3= 0.0; %wavelength/2; %input('Enter position error in meters 

: '); 

fail3=0.0; %5.0; %input('Enter number of elements failed (must be 

less than MxN) :'); 

theta0=0.0; %pi*rand(1,1)-pi/2; %input('Enter of AOA elevation angle 

in radians (-pi/2 < angle < pi/2) :'); 
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thetal=0.0; %theta0-pi/10 % pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) 

theta2=0.0; %theta0+pi/10 %pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) :'); 

theta3=0.0; 

phi01=pi/4; %disp('Assume AOA azimuth angle = pi/4 radians (45 degrees) 

'); 

phi02=pi/4; 
phi03=pi/4; 

phil=pi/2; %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

phi2=pi/2; %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 
phi3=pi/2; 

fsl=50; %input('Enter sampling rate (>300) :'); 

fs2=100; 

fs3=500; 


AF_mag_sum_thetal=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta2=zeros(1,3142); 
affect the AF_cal function 

"O 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta3=zeros(1,3142) ; 
affect the AF_cal function 

"O 

Do 

not 

change 

this. 

it 

will 


AF_mag_sum_phil=zeros(1,6284); 

AF_mag_sum_phi2=zeros(1, 6284); 

AF_mag_sum_phi3=zeros(1, 6284) ; 

%load('Array'); 

[x,x_err,y,y_err,distance,midpointX,midpointY]=Random_LOC_DOA(M, N, L, err 
or2) ; 

%x=x_err; 

%y=y_err; 

figure(l), plot(x,y, ' * ' ), grid on, hold on 
for j=l:1:Na 

j 

%[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 

%x (1, 1) =0.0; 

%x_err(1,1)=0.0; 

%y(1,1)=0.0; 

%y_err(1,1)=0.0; 


[mag_Wtsl,ang_Wtsl]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO1,phil,phi2,fs 
1,wavelength,f); 

[mag_Wts_faill]=Random_FAIL(N,M,mag_Wtsl,faill); 

[THETAl, 

AF_mag_thetal]=AF_cal_theta(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
phiO 1) ; 

[PHIl, 

AF_mag_phil]=AF_cal_phi(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
thetaO); 

AF_mag_sum_thetal=AF_mag_sum_thetal+AF_mag_thetal; 

AF_mag_sum_phil=AF_mag_sum_phi1+AF_mag_phi1; 

%end 

avg_AF_mag_t het a 1 =AF_mag_sum_t he t a 1; 
avg_AF_mag_phil=AF_mag_sum_phi1; 

%for j=l:1:Na 
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%[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 

%x (1, 1) =0.0; 

%x_err(1,1)=0.0; 

%y(1,1)=0.0; 

%y_err(1,1)=0.0; 

[mag_Wts2,ang_Wts2]=LMS(N,M,x,y,thetaO,thetal,theta2,phi02,phil,phi2,fs 

2, wavelength,f); 

[mag_Wts_fail2]=Random_FAIL(N,M,mag_Wts2,fail2); 

[THETA2, 

AF_mag_theta2]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,wave 
length, phi02); 

[PHI2, 

AF_mag_phi2]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,waveleng 
th, thetaO); 

AF_mag_sum_theta2=AF_mag_sum_theta2+AF_mag_theta2; 
AF_mag_sum_phi2=AF_mag_sum_phi2+AF_mag_phi2; 

%end 

avg_AF_mag_theta2=AF_mag_sum_theta2; 
avg_AF_mag_phi2=AF_mag_sum_phi2; 

%for j=l:1:Na 

%[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 

%x (1, 1) =0.0; 

%x_err(1,1)=0.0; 

%y(1,1)=0.0; 

%y_err(1,1)=0.0; 

[mag_Wts3,ang_Wts3]=LMS(N,M,x,y,thetaO,thetal,theta2,phi03,phil,phi2,fs 

3, wavelength,f); 

[mag_Wts_fail3]=Random_FAIL(N,M,mag_Wts3,fail3); 

[THETA3, 

AF_mag_theta3]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,wave 
length, phi03); 

[PHI3, 

AF_mag_phi3]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,waveleng 
th, thetaO); 

AF_mag_sum_theta3=AF_mag_sum_theta3+AF_mag_theta3; 
AF_mag_sum_phi3=AF_mag_sum_phi3+AF_mag_phi3; 

end 

avg_AF_mag_theta3=AF_mag_sum_theta3; 
avg_AF_mag_phi3=AF_mag_sum_phi3; 
thetal=(THETAl/pi)*180; 
theta2=(THETA2/pi)*180; 
theta3=(THETA3/pi)*180; 
phil=(PHIl/pi)*180; 
phi2=(PHI2/pi)*180; 
phi3=(PHI3/pi)*180; 

% 2D plot 
"6 

figure (2), %subplot (3,1,1) 

plot(thetal,20*logl0(avg_AF_mag_thetal),'r'), hold on, grid on 

%xlabel('Elevation Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), 

%xlim([-40 40]),ylim([-20 25]) 
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%title('50 samples') 

%subplot (3,1,2) 

plot(theta2,20*logl0(avg_AF_mag_theta2),'b'), grid on 

%xlabel('Elevation Angle-of-Arrival (degrees)'), ylabel('Average power 
gain (dB)'), 

%xlim([-40 40]),ylim([-20 25]) 

%title('100 samples'); 

%subplot(3,1,3) 

plot(theta3,20*logl0(avg_AF_mag_theta3),'g'), grid on 

xlabel('Elevation Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), xlim([-40 40]),ylim([0 40]) 

title ('25 Element Sensor Field (20 iterations) ') 

legend('50 samples','100 samples', '500 samples') 

plotedit on 

zoom on 

orient tall 

% 2D Polar plot 


%figure (2), polardbl(THETAl, 20*logl0(avg_AF_mag_thetal),-60, 'm'), grid 
on 

%orient tall 

figure (3), polar(THETAl,20*logl0(avg_AF_mag_thetal),'r') , grid on, 
hold on 

polar(THETA2,20*logl0(avg_AF_mag_theta2),'b'), 
polar(THETA3,20*logl0(avg_AF_mag_theta3),'g') 
axis auto 

legend('50 samples','100 samples', '500 samples') 
title('Adaptive Beamforming - 157 MHz') 

gtext('5\lambda square meter grid, where \lambda = 1.9108 m') 
orient tall% - 


% Filename: Node location errors 



% Input parameters: 

% Na = number of random arrays to be generated 
% M = number of elements along x-axis (M must be equal to N) 

% N = number of elements along y-axis (N must be equal to M) 

% wavelength = design wavelength of array (in meters) 

% L = length of square planar (in meters) 

% error = position errors in meters (assumed uniformly distributed 
between [-posterror and iposterror]) 

% fail = number of elements failure 

% thetaO = AOA theta angle (assume AOA phi angle = pi/4 or 90 degrees) 

% fs = sampling rate for the LMS algorithm 


% Output: 

% AF_mag = array factor 


%function[]=Random_MAIN() 
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clear 


Na=10; %input('Enter number of random array samples to be generated 

M=5; %input('Enter number of elements along x-axis, M (must be >1) 

: '); 

N=M; 
c=3e8; 
f=157e6; 

%disp('Number of elements along y-axis, N = M') 

wavelength=c/f; %input('Enter design wavelength of array in meters 

: '); 

L=50*wavelength; %input('Enter length of sensor cluster square 

area in meters : ' ) ; 

errorl= 0.0; %input('Enter position error in meters 

faill=0.0; %input('Enter number of elements failed (must be less 

than MxN) :'); 

error2= wavelength/2; %input('Enter position error in meters :'); 
fail2=0.0; %input('Enter number of elements failed (must be less 

than MxN) : ') ; 

error3= wavelength/2; %input('Enter position error in meters :'); 
fail3=5.0; %input('Enter number of elements failed (must be less 

than MxN) : ' ) ; 

theta0=pi/4; %pi*rand(1,1)-pi/2; %input('Enter of AOA elevation angle 

in radians (-pi/2 < angle < pi/2) :'); 

thetal=pi/4; %theta0-pi/10 % pi*rand(1,1)-pi/2; %input('Enter of AOA 

elevation angle in radians (-pi/2 < angle < pi/2) :'); 

theta2=pi/4; %theta0+pi/10 %pi*rand(1,1)-pi/2; %input('Enter of AOA 
elevation angle in radians (-pi/2 < angle < pi/2) :'); 

phi01=pi/4 %disp('Assume AOA azimuth angle = pi/4 radians (45 degrees) 

'); 

phi02=pi/4 
phi03=pi/4 

phil=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

phi2=pi/10 %2*pi*rand(1,1)-pi; %disp('Assume AOA azimuth angle = pi/4 
radians (45 degrees) '); 

fs=500; %input('Enter sampling rate (>300) :'); 


AF_mag_sum_thetal=zeros(1,3142) ; 
affect the AF_cal function 

"o 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta2=zeros(1,3142); 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 

AF_mag_sum_theta3=zeros(1,3142) ; 
affect the AF_cal function 

"6 

Do 

not 

change 

this. 

it 

will 


AF_mag_sum_phil=zeros(1,6284); 

AF_mag_sum_phi2=zeros(1, 6284) ; 

AF_mag_sum_phi3=zeros(1, 6284) ; 

%load('Array') ; 

%[x, x_err, y,y_err,distance,midpointX,midpointY]=Random_LOC_DOA(M,N,L,er 
ror2) ; 

% figure(l), plot(x,y, ' * ' ), grid on, hold on 
for j=l:1:Na 

j 

[x,x_err,y,y_err]=Random_LOC(M,N,L,error2); 
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grid on, hold on 


x(l,1)=0.0; 
x_err(1,1)=0.0; 
y(1,1)=0.0; 
y_err(1,1)=0.0; 
figure(1), plot(x,y, 

[mag_Wtsl,ang_Wtsl]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO1, phil, phi2, fs 
,wavelength,f); 

[mag_Wts_faill]=Random_FAIL(N,M,mag_Wtsl, faill); 

[THETAl, 

AF_mag_thetal]=AF_cal_theta(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
phiO 1) ; 

[PHIl, 

AF_mag_phil]=AF_cal_phi(N,M,x,y,mag_Wts_faill,ang_Wtsl,wavelength, 
thetaO); 

AF_mag_sum_t he t a 1 =AF_mag_sum_t het a 1 +AF_mag_t het a 1; 

AF_mag_sum_phil=AF_mag_sum_phi1+AF_mag_phi1; 


[mag_Wts2,ang_Wts2]=LMS(N,M,x,y,thetaO,thetal,theta2,phi02,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail2]=Random_FAIL(N,M,mag_Wts2, fail2); 

[THETA2, 

AF_mag_theta2]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,wave 
length, phi02); 

[PHI2, 

AF_mag_phi2]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail2,ang_Wts2,waveleng 
th, thetaO); 

AF_mag_sum_theta2=AF_mag_sum_theta2+AF_mag_theta2; 
AF_mag_sum_phi2=AF_mag_sum_phi2+AF_mag_phi2; 


[mag_Wts3,ang_Wts3]=LMS(N,M,x,y,thetaO,thetal,theta2,phiO3,phil,phi2,fs 
,wavelength,f); 

[mag_Wts_fail3]=Random_FAIL(N,M,mag_Wts3,fail3); 

[THETA3, 

AF_mag_theta3]=AF_cal_theta(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,wave 
length, phi03); 

[PHI3, 

AF_mag_phi3]=AF_cal_phi(N,M,x_err,y_err,mag_Wts_fail3,ang_Wts3,waveleng 
th, thetaO); 

AF_mag_sum_theta3=AF_mag_sum_theta3+AF_mag_theta3; 
AF_mag_sum_phi3=AF_mag_sum_phi3+AF_mag_phi3; 


end 

xlabel('x-direction in meters (m)'), ylabel('y-direction in meters 
(m) ' ) 

title ('SOMambda square meter grid') 
avg_AF_mag_t het a1=AF_mag_sum_t het a1; 
avg_AF_mag_phil=AF_mag_sum_phi1; 

%avg_AF_mag_theta2=AF_mag_sum_theta2; 
avg_AF_mag_phi2=AF_mag_sum_phi2; 

%avg_AF_mag_theta3=AF_mag_sum_theta3; 
avg_AF_mag_phi3=AF_mag_sum_phi3; 
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%fname = input('Save average array factor data for elevation angle as 
(with .mat extension) : 's'); 

%save(fname,'THETA','avg_AF_mag_theta') 

%fname = input('Save average array factor data for azimuth as (with 
.mat extension) : ', 's'); 

%save(fname,'PHI','avg_AF_mag_phi') 

theta=(THETAl/pi)*180; 
phil=(PHIl/pi)*180; 
phi2=(PHI2/pi)*180; 
phi3=(PHI3/pi)*180; 

theta_degrees=(thetaO/pi)*180 
phi_degreesl=(phiOl/pi)*180 
phi_degrees2=(phi02/pi)*180 
phi_degrees3=(phi03/pi)*180 

% 2D plot 


figure (2), subplot (3,1,1) 

plot(phil, 20*logl0(avg_AF_mag_phil), 'r'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB)'), xlim([0 360]),ylim([-10 30]) 
subplot (3,1,2) 

plot(phi2,20*logl0(avg_AF_mag_phi2),'b'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB) ') , xlim([0 360]),ylim([-10 30]) 
title('l meter error'); 
subplot(3,1,3) 

plot(phi3,20*logl0(avg_AF_mag_phi3),'g'), grid on 

xlabel('Azimuth Angle-of-Arrival (degrees)'), ylabel('Average power 

gain (dB) ') , xlim([0 360]),ylim([-10 30]) 

title('l meter error and 20 nodes fail') 

plotedit on 

zoom on 

orient tall 

% 2D Polar plot 


figure (3), polardbl(THETAl, 20*logl0(avg_AF_mag_thetal),-60, 'm'), grid 
on 

orient tall 

figure (4), polar(PHIl,20*logl0(avg_AF_mag_phil),'r'), grid on, hold on 
polar(PHI2,20*logl0(avg_AF_mag_phi2), 'b') , 
polar(PHI3, 20*logl0(avg_AF_mag_phi3), 'g') 

legend('No errors and no failures','\lambda/2 error in node location', 
'\lambda/2 error in node location and 5 node failures') 
title('Adaptive Beamforming - 157 MHz') 

gtext('50\lambda square meter grid, where \lambda = 1.9108 m') 
orient tall 

%PG_plot(THETA,avg_AF_mag_theta) 

%PG_plot(PHI,avg_AF_mag_phi) 

%PG_plot_polar(THETA,avg_AF_mag) 

%load('AF_256_0deg_avg.mat','THETA','avg_AF_mag'); 
polardbl(THETA,20*logl0(avg_AF_mag),-60, 'm'); 
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APPENDIX Q. CRAMER-RAO BOUND 


Using the example taken from [71] the Cramer-Rao bound will be derived. 


h(x) = £'[x I x] -x; eonditional bias 
First rule; 


(Q-1) 


b{x) = E^x{y)-x\\x]= ^[x{y)-x~\p{y\x)dy 

-00 

00 

cib(x) dE{[yy)-i\u] 
dx dx dx 

^ x[y)p[y\x)d^--^^"^ xp[y\x)d^ 


dx 


Using Leibniz’s Rule: 
db[x) d 


= j ^[_Ky)p{y\^)\^y- j ^[_w{y\x)\dy 


dx dx 

dp{y\x) 


\ Hy f^^2''^ dy+ \ ^^p{y\x)dy- j j ^p{y\x)dy 


dx 


dx 


dx 


dx 
dx ' 


(Q-2) 


db(^x) 

dx 


-J %piy\^)^y^ j ^iy) 


dp[y\x) 


dx 


dy 


-00 -oo 

Chain Rule: 


l /p{y\^) 

■L dx 


dy 


du d\nu , X , . , du , , .d\np{y\x) 

— = u -: where u = p(y\x) and — = P\y\ x - ^ - - 
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-00 -00 

j p{y\x)dy^ j [^(j)-^]^— 

-00 -00 


=1 


j[x(j)-x] 


d\np[y\x) 

dx 


[p{y\x)\dy = 


Shwarz inequality: 


j f{t)g{t)dt 


< 



/(j) = V^(j|x)[x(y)-x] 


After substitution equation (Q-3) into equation (Q-4): 


^^db{x) 

< 

dx 

A 


x)[x(y)-x] Jxj yjp(j\^ 


Jln;7(y |x) 


dx 


< 


i 


00 2 oo 

j|[x(y)-x]| p(ylx)dxj 


dlnp(ylx) 


dx 


p[y\x)dx 


< 


i 


j|[x(y)-x]| p(y\x)dxj 

-00 '' -V--00 

second moment given x 


Jln;7(y |x) 


dx 


p[y\x)dx 


second moment given x 


dx 


Rewriting: 


(Q-3) 


, (Q-4) 


(Q-5) 
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dx 



1 



E 

^ d\np[y\x)^ 

2 T 

X 

Jj 


^ dx j 



{x{y)-xf 



1 I 

dx 


77 

^ d\np[y\x)^ 

2 T 

J2j 

^ dx j 

1 X 


E 




^^db{x) 

dx 

2 

p 

( d\np[y\xy 

2 

V 


1 

dx 

7 

1 Jk 


(Q-6) 


For an unconditional bias estsimator, which is used in this research 

b = E[x-x] = Q 

And for a known density the eondition on x is removed. Then the above reduees to the 
following; 


(x(y)-x)' 


> 


1 


d\np[y\x) 

dx 


As per Whalen, the equality in the Schwarz inequality holds only when 

f[t) = kg{t) 


(Q-V) 


(x(y)-x)' 


1 


d\np[y\x) 

dx 




( d\n pi^y\xy^ 

1 X 

_ rf dp{y\xT 

y dx j 

J 

1 dx , 

-<x> V y 


p ' {y I x)dy 


(Q-8) 


267 



Continuing with example from Whalen, using chain rule equation (Q-8) can be simplified 
to the familiar form of equation (Q-9) 


u=p{y\x)-, 

d\nu _ u^^du 
dx dx 


du _ud\nu 
dx dx 


E — \np{y\x) =\ 


j P{y\x)dy- j 


d^p{y\x) (dp{y\x) 


dp{y\x) 


P ' {y I x) \dy 


P^'{y\x)dy 


d\np[y\x) 


E\\x[y)-x~\ 


X =-E 


d\np[y\x) 


d^\np[y\x) 


(Q-9) 


Now to use this information to estimate the Cramer-Rao bound to estmate an unknown 
parameter. 


r[t) = s^[t\x) + n[t)-. 


^0 - ^ - ^0 + ^ ’ 


(Q-10) 


It is assumed that the signal is known over the same time interval. The probability density 
function of the received signal in white Gaussian noise is as follows: 


^samp 1 1 

"TT 1 


p{r.)= n 

i=0 ■sJZTTO'. 


N =■ 

samp 


sampling interval 


(Q-11) 


2 
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p{l) = 


\yl2^a j 


,_n 2 (J 


p{l) = 


( 1 A 


V J 

N 

samp 


N -1 
‘'samp ‘ 2 


■=0 


' T = 

1 


23, 

2 2 




v^^oy 


=0 2 - 


"■'“P '(»;.-f;.f<Si 


(Q-12) 


r{t) = E[r{t)\- 
:. p{r^ = ^exp 
( N 

K(rh ^ Hr) 

\ ^ J 


Let 5 ^ 00 , so that dt^Q 


] (n-nf 
J M 

0 ‘'o 


dt 


(Q-13) 


E[r{t)\ = Sr{t 

II); 


- V - 

II 

1^4 exp 

J N 

0 '^'^O J 


In j!7(r I x) = lny4 + In 


exp 


silEf 

J M 

0 


dt 


d 

L L 

r (r. — 

‘ '> dt 

51n/>(r x) 

0 ’0 


dx, 


5x 
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Substitute equation (Q-10) into equation (Q-13): 


51n/»(r I x) 


I 

-i 




0 ^^0 


dt 


dx^ 

Na dx, 


5x. 

T 

\(r{t)-s^{t\x)f dt 


; substitute in for the reeeived signals 


(Q-14) 


Now taking the partial derivatives 


j2(r(t)-5^ [t I x))' [t I x)\dt 

5x. A^oo 


r{t) = s^{t\x)-n{t) 

|x) = n(t); where 



The noiseless eomplex envelop of the reeeived signal: 

{t) = Aa{t)cosi^o}J + (l){t) + 6^^ 
r (t) = (t)exp(76*^) + n(t); where, 


(Q-15) 


(Q-16) 


As per Whalen, 6c is eonsidered a nuisanee parameter of Uniform Density and will be 
removed by taking the expeetation given 6c. 


^P(0l^c]= j ^exp 

71 

= J Aexp 


jfzkWl,, 




iVn 


iVn 


1 


{6)d6; where {6) = ^ 

ITT 


-^6 

2n 


(Q-17) 


In 


71 

\ 


exp 




iVn 


d6 
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Now substitute equation (Q-16) into equation (Q-17) as follows: 




In 


0 ^^0 


d0\ where 


1 ^ 

5 (t) = y 4 a(t)exp| 7 [^zi(t) + ^^]| and —^a^{t^dt = \ 





T 

f y 


Jexp 

--Lj 

N •* 

r^^-2r[t)s^{t) + s^j^ 

dt 

-n 

”0 0 

vho^'f) *fh'(oy 



de 


A 

2n 


JL 

jexp 


|f(t)|"-2^a(t)f(t)exp{7[^(t) + ^]}+ |^|" |a(t)p 

^0 ol L =1 jy 


dt 


de 


A_ 

2n 


exp 


^ \r{t^ +Uf \a{t^'‘^ 


N, 


0 0 


dt 


J — J^2y4a(t)f (t)exp|7[^zi(t) + ^]|j 

-;r L^O 0 

, (Q-18) 

For eonvienee Whalen redefines the last integral term as follows: 


de\ 


1 

J a {t)?{t) exp [ 7 ^(t)] dt = 2q exp ( 7 ^ 0 ) 


A 


E[r{t)\e~] = —\^^V 


-j-]\^?{tf +\Af\a{t)\]^dt Jexp 

-''^00 -71 


I [ 2 ^^ exp ( 7 ^) exp ( 7 ^)] ^^ 




de 


, (Q-19) 


Where the modified Bessel funetion is defined as follows: 


2 jexp [/eos ( 6 *)] <76* = 2 ;r/o (f) 




^ '^\Htf+\A[\a{e'" 


N, 


0 0 


dt 


exp 


I [2Aq exp {j/3 ) exp ( 76 >)] 




de] 


(Q-20) 
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E[f(l)lff] = 2^£/„ 


^2Aq^ 

V^y 


exp 


JV, 


0 0 


= (^)|'+|^|^|a(^)f 


dt 


dropping the given 6 term 


E\r{t)yAI, 


^2Aq^ 

V^o y 


exp 


0 0 


-11 VitW + \A\ \a(t 

N, 


dt 


In^ f(t) =ln<y4/o 


^2Aq'^ 

V^o y 


exp 


N, 


0 0 


rltll +1^41^ aid' 


dt 


I |2 

remembering that a(t) = (t) 


In^ r{t) =lny4 + ln 


^ 2 Aq'^ 
V^oy 


+ In I exp 


N, 


0 0 


^(or+MrK(o 


dt 


for high SNR: f (t) = Ad^ defined earlier less the noise term 

for a large argument the natural log of the modified Bessel funetion is as follows: 

In/g (x) = X 


In^ =lnA + 


2Aq 


+ 


(0< +\Af\d id'^ 

0 


dt 


1 A 

= In ^ +- - + 






= ln^ + 


2Aq 2A 


0 0 
2 T 


[t 


dt 


No N, , 


[t 


dt 


(Q-21) 


, (Q-22) 
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